{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ac4e1c41",
   "metadata": {},
   "source": [
    "# 13 - 面板数据和固定效应\n",
    "\n",
    "## 控制你看不到的东西\n",
    "\n",
    "倾向得分、线性回归和匹配等方法非常擅长控制非随机数据中的混淆现象，但它们依赖于一个关键假设：控制后无混淆\n",
    "\n",
    "$\n",
    "(Y_0, Y_1) \\perp T | X\n",
    "$\n",
    "\n",
    "简而言之，它们要求所有混淆因素都是已知的和可测量的，这样我们才能以它们为条件并使处理尽可能随机。一个主要问题是有时我们根本无法衡量混淆因素。例如，拿一个经典的劳动经济学问题来计算婚姻对男性收入的影响。经济学中众所周知的事实是，已婚男性的收入高于单身男性。但是，尚不清楚这种关系是否是因果关系。可能是受过更多教育的男性更有可能结婚，也更有可能从事高收入工作，这意味着教育是婚姻对收入影响的一个混淆因素。对于这个混淆因素，我们可以衡量研究中人员的教育程度，并对其进行回归控制。但另一个混淆因素可能是好看的外貌。可能是更英俊的男人更有可能结婚，也更有可能获得高薪工作。不幸的是，美丽和智力这类特征类似，这是我们无法很好衡量的东西。\n",
    "\n",
    "这使我们陷入困境，因为如果我们有无法衡量的混淆因素，我们就有偏差。正如我们之前所见，处理这个问题的一种方法是使用工具变量。但是想出好的工具变量并不是一件容易的事，需要很多的创造力。在这里，让我们看一下利用时间或数据的时间结构的替代方案。\n",
    "\n",
    "这个想法是使用**面板数据**。当我们**在多个时间段内对同一个人进行观察**，我们就得到了面板数据。面板数据格式在行业中非常常见，比如当对同一用户在多个时间段的行为记录进行记录和保存的时候。我们可以利用面板数据的原因是因为我们可以在干预前后比较同一个单位，看看他们在干预前后的表现如何。在我们深入研究数学之前，让我们看看如何进行直观地理解。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "863afdfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import style\n",
    "from matplotlib import pyplot as plt\n",
    "import statsmodels.formula.api as smf\n",
    "import graphviz as gr\n",
    "\n",
    "from linearmodels.panel import PanelOLS\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "pd.set_option(\"display.max_columns\", 6)\n",
    "style.use(\"fivethirtyeight\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8fc10b2",
   "metadata": {},
   "source": [
    "首先，让我们看一下我们在跨时间包含同一单位的多个观察后得到的因果图。 假设我们有这样一种情况，第一次结婚同时导致收入和随后的婚姻状况。 对于第 2 次和第 3 次也是如此。此外，假设所有时间段的美貌程度都是相同的（这是一个大胆的说法，但如果时间只有几年，这是合理的），它会导致婚姻和收入。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "e220106a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"255pt\" height=\"332pt\"\r\n",
       " viewBox=\"0.00 0.00 254.79 332.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 328)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-328 250.794,-328 250.794,4 -4,4\"/>\r\n",
       "<!-- 婚姻_1 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>婚姻_1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"102.397\" cy=\"-234\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"102.397\" y=\"-230.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">婚姻_1</text>\r\n",
       "</g>\r\n",
       "<!-- 收入_1 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>收入_1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"49.3968\" cy=\"-162\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"49.3968\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">收入_1</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_1&#45;&gt;收入_1 -->\r\n",
       "<g id=\"edge1\" class=\"edge\"><title>婚姻_1&#45;&gt;收入_1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M90.1043,-216.765C83.4819,-208.018 75.1789,-197.052 67.8017,-187.308\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"70.4494,-185.007 61.6226,-179.147 64.8686,-189.233 70.4494,-185.007\"/>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_2 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>婚姻_2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"140.397\" cy=\"-162\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"140.397\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">婚姻_2</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_1&#45;&gt;婚姻_2 -->\r\n",
       "<g id=\"edge2\" class=\"edge\"><title>婚姻_1&#45;&gt;婚姻_2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M111.402,-216.411C115.945,-208.042 121.554,-197.71 126.625,-188.37\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"129.761,-189.928 131.456,-179.47 123.609,-186.589 129.761,-189.928\"/>\r\n",
       "</g>\r\n",
       "<!-- 收入_2 -->\r\n",
       "<g id=\"node4\" class=\"node\"><title>收入_2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"36.3968\" cy=\"-90\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"36.3968\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">收入_2</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_2&#45;&gt;收入_2 -->\r\n",
       "<g id=\"edge3\" class=\"edge\"><title>婚姻_2&#45;&gt;收入_2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M119.351,-146.834C103.885,-136.424 82.6742,-122.148 65.5728,-110.638\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"67.3067,-107.586 57.0565,-104.906 63.3981,-113.393 67.3067,-107.586\"/>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_3 -->\r\n",
       "<g id=\"node5\" class=\"node\"><title>婚姻_3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.397\" cy=\"-90\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"178.397\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">婚姻_3</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_2&#45;&gt;婚姻_3 -->\r\n",
       "<g id=\"edge4\" class=\"edge\"><title>婚姻_2&#45;&gt;婚姻_3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M149.402,-144.411C153.945,-136.042 159.554,-125.71 164.625,-116.37\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"167.761,-117.928 169.456,-107.47 161.609,-114.589 167.761,-117.928\"/>\r\n",
       "</g>\r\n",
       "<!-- 收入_3 -->\r\n",
       "<g id=\"node6\" class=\"node\"><title>收入_3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"210.397\" cy=\"-18\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"210.397\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">收入_3</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_3&#45;&gt;收入_3 -->\r\n",
       "<g id=\"edge5\" class=\"edge\"><title>婚姻_3&#45;&gt;收入_3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M186.143,-72.055C189.838,-63.9726 194.341,-54.1214 198.46,-45.1117\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"201.747,-46.3389 202.722,-35.789 195.381,-43.4286 201.747,-46.3389\"/>\r\n",
       "</g>\r\n",
       "<!-- 美貌程度 -->\r\n",
       "<g id=\"node7\" class=\"node\"><title>美貌程度</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"134.397\" cy=\"-306\" rx=\"44.393\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"134.397\" y=\"-302.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">美貌程度</text>\r\n",
       "</g>\r\n",
       "<!-- 美貌程度&#45;&gt;婚姻_1 -->\r\n",
       "<g id=\"edge6\" class=\"edge\"><title>美貌程度&#45;&gt;婚姻_1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M126.65,-288.055C122.956,-279.973 118.452,-270.121 114.334,-261.112\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"117.413,-259.429 110.072,-251.789 111.046,-262.339 117.413,-259.429\"/>\r\n",
       "</g>\r\n",
       "<!-- 美貌程度&#45;&gt;收入_1 -->\r\n",
       "<g id=\"edge9\" class=\"edge\"><title>美貌程度&#45;&gt;收入_1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M104.921,-292.561C88.0338,-283.804 68.1222,-270.356 57.3968,-252 46.5361,-233.413 45.099,-208.826 46.0604,-190.205\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"49.5664,-190.262 46.8505,-180.021 42.5873,-189.72 49.5664,-190.262\"/>\r\n",
       "</g>\r\n",
       "<!-- 美貌程度&#45;&gt;婚姻_2 -->\r\n",
       "<g id=\"edge7\" class=\"edge\"><title>美貌程度&#45;&gt;婚姻_2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M139.671,-288.007C142.59,-277.704 145.919,-264.227 147.397,-252 149.891,-231.357 147.849,-207.908 145.34,-190.263\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"148.754,-189.444 143.754,-180.104 141.838,-190.524 148.754,-189.444\"/>\r\n",
       "</g>\r\n",
       "<!-- 美貌程度&#45;&gt;收入_2 -->\r\n",
       "<g id=\"edge10\" class=\"edge\"><title>美貌程度&#45;&gt;收入_2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M104.14,-292.548C84.9439,-283.491 60.5771,-269.721 43.3968,-252 18.0646,-225.871 13.3374,-215.278 4.39677,-180 0.466042,-164.49 0.104481,-159.414 4.39677,-144 7.17887,-134.01 12.4094,-124.021 17.8476,-115.453\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"20.7954,-117.341 23.496,-107.096 14.9961,-113.421 20.7954,-117.341\"/>\r\n",
       "</g>\r\n",
       "<!-- 美貌程度&#45;&gt;婚姻_3 -->\r\n",
       "<g id=\"edge8\" class=\"edge\"><title>美貌程度&#45;&gt;婚姻_3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M144.144,-288.346C157.074,-265.054 179.009,-220.991 185.397,-180 188.599,-159.455 186.522,-136.001 183.82,-118.33\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"187.218,-117.428 182.099,-108.152 180.316,-118.595 187.218,-117.428\"/>\r\n",
       "</g>\r\n",
       "<!-- 美貌程度&#45;&gt;收入_3 -->\r\n",
       "<g id=\"edge11\" class=\"edge\"><title>美貌程度&#45;&gt;收入_3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M152.558,-289.462C162.952,-279.674 175.464,-266.229 183.397,-252 215.74,-193.983 214.426,-173.815 223.397,-108 225.558,-92.1466 225.316,-87.8844 223.397,-72 222.358,-63.4026 220.403,-54.1878 218.32,-45.8926\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"221.637,-44.748 215.671,-35.9927 214.875,-46.5576 221.637,-44.748\"/>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1f97c6576c8>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = gr.Digraph()\n",
    "g.edge(\"婚姻_1\", \"收入_1\")\n",
    "g.edge(\"婚姻_1\", \"婚姻_2\")\n",
    "g.edge(\"婚姻_2\", \"收入_2\")\n",
    "g.edge(\"婚姻_2\", \"婚姻_3\")\n",
    "g.edge(\"婚姻_3\", \"收入_3\")\n",
    "\n",
    "g.edge(\"美貌程度\", \"婚姻_1\")\n",
    "g.edge(\"美貌程度\", \"婚姻_2\")\n",
    "g.edge(\"美貌程度\", \"婚姻_3\")\n",
    "\n",
    "g.edge(\"美貌程度\", \"收入_1\")\n",
    "g.edge(\"美貌程度\", \"收入_2\")\n",
    "g.edge(\"美貌程度\", \"收入_3\")\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7df33444",
   "metadata": {},
   "source": [
    "请记住，我们无法控制外貌多美这个因素，因为我们无法衡量它。 但是我们仍然可以使用面板结构，所以它不再是问题了。 这个想法是，我们可以将美——以及任何其他随时间不变的属性——视为定义一个人的综合要素。 虽然我们不能直接控制它们，但我们可以控制个人本身。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a9e30326",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
       "<!-- Generated by graphviz version 2.38.0 (20140413.2041)\r\n",
       " -->\r\n",
       "<!-- Title: %3 Pages: 1 -->\r\n",
       "<svg width=\"255pt\" height=\"332pt\"\r\n",
       " viewBox=\"0.00 0.00 254.79 332.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 328)\">\r\n",
       "<title>%3</title>\r\n",
       "<polygon fill=\"white\" stroke=\"none\" points=\"-4,4 -4,-328 250.794,-328 250.794,4 -4,4\"/>\r\n",
       "<!-- 婚姻_1 -->\r\n",
       "<g id=\"node1\" class=\"node\"><title>婚姻_1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"102.397\" cy=\"-234\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"102.397\" y=\"-230.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">婚姻_1</text>\r\n",
       "</g>\r\n",
       "<!-- 收入_1 -->\r\n",
       "<g id=\"node2\" class=\"node\"><title>收入_1</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"49.3968\" cy=\"-162\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"49.3968\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">收入_1</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_1&#45;&gt;收入_1 -->\r\n",
       "<g id=\"edge1\" class=\"edge\"><title>婚姻_1&#45;&gt;收入_1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M90.1043,-216.765C83.4819,-208.018 75.1789,-197.052 67.8017,-187.308\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"70.4494,-185.007 61.6226,-179.147 64.8686,-189.233 70.4494,-185.007\"/>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_2 -->\r\n",
       "<g id=\"node3\" class=\"node\"><title>婚姻_2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"140.397\" cy=\"-162\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"140.397\" y=\"-158.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">婚姻_2</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_1&#45;&gt;婚姻_2 -->\r\n",
       "<g id=\"edge2\" class=\"edge\"><title>婚姻_1&#45;&gt;婚姻_2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M111.402,-216.411C115.945,-208.042 121.554,-197.71 126.625,-188.37\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"129.761,-189.928 131.456,-179.47 123.609,-186.589 129.761,-189.928\"/>\r\n",
       "</g>\r\n",
       "<!-- 收入_2 -->\r\n",
       "<g id=\"node4\" class=\"node\"><title>收入_2</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"36.3968\" cy=\"-90\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"36.3968\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">收入_2</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_2&#45;&gt;收入_2 -->\r\n",
       "<g id=\"edge3\" class=\"edge\"><title>婚姻_2&#45;&gt;收入_2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M119.351,-146.834C103.885,-136.424 82.6742,-122.148 65.5728,-110.638\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"67.3067,-107.586 57.0565,-104.906 63.3981,-113.393 67.3067,-107.586\"/>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_3 -->\r\n",
       "<g id=\"node5\" class=\"node\"><title>婚姻_3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"178.397\" cy=\"-90\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"178.397\" y=\"-86.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">婚姻_3</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_2&#45;&gt;婚姻_3 -->\r\n",
       "<g id=\"edge4\" class=\"edge\"><title>婚姻_2&#45;&gt;婚姻_3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M149.402,-144.411C153.945,-136.042 159.554,-125.71 164.625,-116.37\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"167.761,-117.928 169.456,-107.47 161.609,-114.589 167.761,-117.928\"/>\r\n",
       "</g>\r\n",
       "<!-- 收入_3 -->\r\n",
       "<g id=\"node6\" class=\"node\"><title>收入_3</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"210.397\" cy=\"-18\" rx=\"36.2938\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"210.397\" y=\"-14.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">收入_3</text>\r\n",
       "</g>\r\n",
       "<!-- 婚姻_3&#45;&gt;收入_3 -->\r\n",
       "<g id=\"edge5\" class=\"edge\"><title>婚姻_3&#45;&gt;收入_3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M186.143,-72.055C189.838,-63.9726 194.341,-54.1214 198.46,-45.1117\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"201.747,-46.3389 202.722,-35.789 195.381,-43.4286 201.747,-46.3389\"/>\r\n",
       "</g>\r\n",
       "<!-- Person (美貌程度, 智力...) -->\r\n",
       "<g id=\"node7\" class=\"node\"><title>Person (美貌程度, 智力...)</title>\r\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"134.397\" cy=\"-306\" rx=\"106.681\" ry=\"18\"/>\r\n",
       "<text text-anchor=\"middle\" x=\"134.397\" y=\"-302.3\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">Person (美貌程度, 智力...)</text>\r\n",
       "</g>\r\n",
       "<!-- Person (美貌程度, 智力...)&#45;&gt;婚姻_1 -->\r\n",
       "<g id=\"edge6\" class=\"edge\"><title>Person (美貌程度, 智力...)&#45;&gt;婚姻_1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M126.487,-287.697C122.79,-279.609 118.311,-269.812 114.223,-260.869\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"117.336,-259.262 109.996,-251.622 110.97,-262.172 117.336,-259.262\"/>\r\n",
       "</g>\r\n",
       "<!-- Person (美貌程度, 智力...)&#45;&gt;收入_1 -->\r\n",
       "<g id=\"edge9\" class=\"edge\"><title>Person (美貌程度, 智力...)&#45;&gt;收入_1</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M98.2208,-288.914C83.0206,-280.231 66.7304,-267.974 57.3968,-252 46.5361,-233.413 45.099,-208.826 46.0604,-190.205\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"49.5664,-190.262 46.8505,-180.021 42.5873,-189.72 49.5664,-190.262\"/>\r\n",
       "</g>\r\n",
       "<!-- Person (美貌程度, 智力...)&#45;&gt;婚姻_2 -->\r\n",
       "<g id=\"edge7\" class=\"edge\"><title>Person (美貌程度, 智力...)&#45;&gt;婚姻_2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M139.671,-288.007C142.59,-277.704 145.919,-264.227 147.397,-252 149.891,-231.357 147.849,-207.908 145.34,-190.263\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"148.754,-189.444 143.754,-180.104 141.838,-190.524 148.754,-189.444\"/>\r\n",
       "</g>\r\n",
       "<!-- Person (美貌程度, 智力...)&#45;&gt;收入_2 -->\r\n",
       "<g id=\"edge10\" class=\"edge\"><title>Person (美貌程度, 智力...)&#45;&gt;收入_2</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M97.1174,-289.131C79.269,-280.171 58.5409,-267.621 43.3968,-252 18.0646,-225.871 13.3374,-215.278 4.39677,-180 0.466042,-164.49 0.104481,-159.414 4.39677,-144 7.17887,-134.01 12.4094,-124.021 17.8476,-115.453\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"20.7954,-117.341 23.496,-107.096 14.9961,-113.421 20.7954,-117.341\"/>\r\n",
       "</g>\r\n",
       "<!-- Person (美貌程度, 智力...)&#45;&gt;婚姻_3 -->\r\n",
       "<g id=\"edge8\" class=\"edge\"><title>Person (美貌程度, 智力...)&#45;&gt;婚姻_3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M144.364,-287.949C157.321,-264.556 179.045,-220.759 185.397,-180 188.599,-159.455 186.522,-136.001 183.82,-118.33\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"187.218,-117.428 182.099,-108.152 180.316,-118.595 187.218,-117.428\"/>\r\n",
       "</g>\r\n",
       "<!-- Person (美貌程度, 智力...)&#45;&gt;收入_3 -->\r\n",
       "<g id=\"edge11\" class=\"edge\"><title>Person (美貌程度, 智力...)&#45;&gt;收入_3</title>\r\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M154.01,-288.084C164.098,-278.431 175.83,-265.572 183.397,-252 215.74,-193.983 214.426,-173.815 223.397,-108 225.558,-92.1466 225.316,-87.8844 223.397,-72 222.358,-63.4026 220.403,-54.1878 218.32,-45.8926\"/>\r\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"221.637,-44.748 215.671,-35.9927 214.875,-46.5576 221.637,-44.748\"/>\r\n",
       "</g>\r\n",
       "</g>\r\n",
       "</svg>\r\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1f97c679308>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = gr.Digraph()\n",
    "g.edge(\"婚姻_1\", \"收入_1\")\n",
    "g.edge(\"婚姻_1\", \"婚姻_2\")\n",
    "g.edge(\"婚姻_2\", \"收入_2\")\n",
    "g.edge(\"婚姻_2\", \"婚姻_3\")\n",
    "g.edge(\"婚姻_3\", \"收入_3\")\n",
    "\n",
    "g.edge(\"Person (美貌程度, 智力...)\", \"婚姻_1\")\n",
    "g.edge(\"Person (美貌程度, 智力...)\", \"婚姻_2\")\n",
    "g.edge(\"Person (美貌程度, 智力...)\", \"婚姻_3\")\n",
    "\n",
    "g.edge(\"Person (美貌程度, 智力...)\", \"收入_1\")\n",
    "g.edge(\"Person (美貌程度, 智力...)\", \"收入_2\")\n",
    "g.edge(\"Person (美貌程度, 智力...)\", \"收入_3\")\n",
    "\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "281057cb",
   "metadata": {},
   "source": [
    "想想看。我们无法衡量美貌和智力等属性，但我们知道拥有这些属性的人是同一个人。因此，我们可以创建一个表示该人的虚拟变量并将其添加到线性模型中。当我们说我们可以控制人本身时，这就是我们的意思：我们正在添加一个表示该特定人的变量（在本例中为虚拟变量）。在我们的模型中使用这个人来估计婚姻对收入的影响时，回归发现婚姻的影响**同时保持人这个变量固定**。添加这个实体个人对应的虚拟变量就是我们所说的固定效应模型。\n",
    "\n",
    "\n",
    "## 固定效应\n",
    "\n",
    "为了方面后面更正式地讲述，让我们首先看一下我们拥有的数据。按照我们的例子，我们将尝试估计婚姻对收入的影响。我们的数据包含多年以来多个个体 (`nr`) 的这两个变量，`married` 和`lwage`。请注意，工资采用对数形式。除此之外，我们还有其他控制措施，例如当年的工作小时数、受教育年限等。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a58dc141",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>nr</th>\n",
       "      <th>year</th>\n",
       "      <th>black</th>\n",
       "      <th>...</th>\n",
       "      <th>lwage</th>\n",
       "      <th>expersq</th>\n",
       "      <th>occupation</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>13</td>\n",
       "      <td>1980</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.197540</td>\n",
       "      <td>1</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>13</td>\n",
       "      <td>1981</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.853060</td>\n",
       "      <td>4</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>13</td>\n",
       "      <td>1982</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.344462</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>13</td>\n",
       "      <td>1983</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.433213</td>\n",
       "      <td>16</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>13</td>\n",
       "      <td>1984</td>\n",
       "      <td>0</td>\n",
       "      <td>...</td>\n",
       "      <td>1.568125</td>\n",
       "      <td>25</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 12 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   nr  year  black  ...     lwage  expersq  occupation\n",
       "0  13  1980      0  ...  1.197540        1           9\n",
       "1  13  1981      0  ...  1.853060        4           9\n",
       "2  13  1982      0  ...  1.344462        9           9\n",
       "3  13  1983      0  ...  1.433213       16           9\n",
       "4  13  1984      0  ...  1.568125       25           5\n",
       "\n",
       "[5 rows x 12 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from linearmodels.datasets import wage_panel\n",
    "data = wage_panel.load()\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba09c474",
   "metadata": {},
   "source": [
    "通常，固定效应模型定义为\n",
    "\n",
    "$\n",
    "y_{it} = \\beta X_{it} + \\gamma U_i + e_{it}\n",
    "$\n",
    "\n",
    "其中 \\\\(y_{it}\\\\) 是个体 \\\\(i\\\\) 在时间 \\\\(t\\\\) 的结果，\\\\(X_{it}\\\\) 是个体变量的向量\\\\(i\\\\) 在时间 \\\\(t\\\\)。 \\\\(U_i\\\\) 是单个 \\\\(i\\\\) 的一组不可观测值。请注意，这些不可观测值随着时间的推移是不变的，因此缺少时间下标。最后，\\\\(e_{it}\\\\) 是错误项。对于教育示例，\\\\(y_{it}\\\\) 是对数工资，\\\\(X_{it}\\\\) 是随时间变化的可观察变量，例如婚姻和经验，\\\\(U_i\\\\)是每个人没有观察到但不变的变量，例如美丽和智力。\n",
    "\n",
    "\n",
    "现在，请记住我说过使用具有固定效果模型的面板数据就像为实体添加虚拟对象一样简单。这是真的，但在实践中，我们实际上并没有这样做。想象一个我们有 100 万客户的数据集。如果我们为它们中的每一个添加一个 dummy，我们最终会得到 100 万列，这可能不是一个好主意。相反，我们使用将线性回归划分为 2 个独立模型的技巧。我们以前见过这个，但现在是回顾它的好时机。假设您有一个线性回归模型，其中包含一组特征 \\\\(X_1\\\\) 和另一组特征 \\\\(X_2\\\\)。\n",
    "\n",
    "$\n",
    "\\hat{Y} = \\hat{\\beta_1} X_1 + \\hat{\\beta_2} X_2\n",
    "$\n",
    "\n",
    "其中 \\\\(X_1\\\\) 和 \\\\(X_1\\\\) 是特征矩阵（每个特征一行，每个观察一列）和 \\\\(\\hat{\\beta_1}\\\\) 和 \\\\(\\hat{ \\beta_2}\\\\) 是行向量。您可以通过执行获得完全相同的 \\\\(\\hat{\\beta_1}\\\\) 参数\n",
    "\n",
    "1. 在第二组特征 \\\\(\\hat{y^*} = \\hat{\\gamma_1} X_2\\\\) 上回归结果 \\\\(y\\\\)\n",
    "2. 在第二个 \\\\(\\hat{X_1} = \\hat{\\gamma_2} X_2\\\\) 上回归第一组特征\n",
    "3. 得到残差 \\\\(\\tilde{X}_1 = X_1 - \\hat{X_1}\\\\) 和 \\\\(\\tilde{y}_1 = y_1 - \\hat{y^*}\\\\)\n",
    "4. 将结果的残差回归到特征残差 \\\\(\\hat{y} = \\hat{\\beta_1} \\tilde{X}_1\\\\)\n",
    "\n",
    "最后一次回归的参数将与使用所有特征运行回归完全相同。但这究竟对我们有什么帮助呢？好吧，我们可以将带有实体假人的模型的估计分解为 2。首先，我们使用假人来预测结果和特征。这些是上面的步骤 1 和 2。\n",
    "\n",
    "现在，还记得在虚拟变量上运行回归是如何像估计该虚拟变量的平均值一样简单吗？如果你不这样做，让我们用我们的数据来证明这是真的。让我们运行一个模型，我们将工资预测为虚拟年份的函数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a092d1eb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>    1.3935</td> <td>    0.022</td> <td>   63.462</td> <td> 0.000</td> <td>    1.350</td> <td>    1.437</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1981]</th> <td>    0.1194</td> <td>    0.031</td> <td>    3.845</td> <td> 0.000</td> <td>    0.059</td> <td>    0.180</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1982]</th> <td>    0.1782</td> <td>    0.031</td> <td>    5.738</td> <td> 0.000</td> <td>    0.117</td> <td>    0.239</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1983]</th> <td>    0.2258</td> <td>    0.031</td> <td>    7.271</td> <td> 0.000</td> <td>    0.165</td> <td>    0.287</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1984]</th> <td>    0.2968</td> <td>    0.031</td> <td>    9.558</td> <td> 0.000</td> <td>    0.236</td> <td>    0.358</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1985]</th> <td>    0.3459</td> <td>    0.031</td> <td>   11.140</td> <td> 0.000</td> <td>    0.285</td> <td>    0.407</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1986]</th> <td>    0.4062</td> <td>    0.031</td> <td>   13.082</td> <td> 0.000</td> <td>    0.345</td> <td>    0.467</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(year)[T.1987]</th> <td>    0.4730</td> <td>    0.031</td> <td>   15.232</td> <td> 0.000</td> <td>    0.412</td> <td>    0.534</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = smf.ols(\"lwage ~ C(year)\", data=data).fit()\n",
    "mod.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33626c7f",
   "metadata": {},
   "source": [
    "请注意该模型如何预测 1980 年的平均收入为 1.3935，1981 年的平均收入为 1.5129 (1.3935+0.1194) 等等。 现在，如果我们按年份计算平均值，我们会得到完全相同的结果。 （请记住，基准年 1980 是截距。因此，您必须将截距添加到其他年份的参数中才能获得该年的平均`lwage`）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fdde9c8f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "year\n",
       "1980    1.393477\n",
       "1981    1.512867\n",
       "1982    1.571667\n",
       "1983    1.619263\n",
       "1984    1.690295\n",
       "1985    1.739410\n",
       "1986    1.799719\n",
       "1987    1.866479\n",
       "Name: lwage, dtype: float64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.groupby(\"year\")[\"lwage\"].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bbaeaf8",
   "metadata": {},
   "source": [
    "这意味着，如果我们得到面板中每个人的平均值，我们基本上是在对其他变量进行个体虚拟回归。这激发了以下估计过程：\n",
    "\n",
    "1. 通过减去个人的平均值来创建时间贬损变量：\n",
    "$\\ddot{Y}_{it} = Y_{it} - \\bar{Y}_i$\n",
    "$\\ddot{X}_{it} = X_{it} - \\bar{X}_i$\n",
    "\n",
    "2. 在 $\\ddot{X}_{it}$ 上回归 $\\ddot{Y}_{it}$\n",
    "\n",
    "\n",
    "请注意，当我们这样做时，未观察到的 \\\\(U_i\\\\) 消失了。由于 \\\\(U_i\\\\) 在时间上是恒定的，所以我们有 \\\\(\\bar{U_i}=U_i\\\\)。如果我们有以下两个方程组\n",
    "\n",
    "$$\n",
    "\\开始{对齐}\n",
    "Y_{it} & = \\beta X_{it} + \\gamma U_i + e_{it} \\\\\n",
    "\\bar{Y}_{i} & = \\beta \\bar{X}_{it} + \\gamma \\bar{U}_i + \\bar{e}_{it} \\\\\n",
    "\\结束{对齐}\n",
    "$$\n",
    "\n",
    "我们从另一个中减去一个，我们得到\n",
    "\n",
    "$$\n",
    "\\开始{对齐}\n",
    "(Y_{it} - \\bar{Y}_{i}) & = (\\beta X_{it} - \\beta \\bar{X}_{it}) + (\\gamma U_i - \\gamma U_i) + ( e_{it}-\\bar{e}_{it}) \\\\\n",
    "(Y_{it} - \\bar{Y}_{i}) & = \\beta(X_{it} - \\bar{X}_{it}) + (e_{it}-\\bar{e}_{它}） \\\\\n",
    "\\ddot{Y}_{it} & = \\beta \\ddot{X}_{it} + \\ddot{e}_{it} \\\\\n",
    "\\结束{对齐}\n",
    "$$\n",
    "\n",
    "它消除了所有未观察到的随时间不变的事物。老实说，不仅未观察到的变量消失了。这发生在所有时间不变的变量上。因此，您不能包含任何随时间保持不变的变量，因为它们将是虚拟变量的线性组合，并且模型不会运行。\n",
    "\n",
    "![img](./data/img/fixed-effects/demeaned.png)\n",
    "\n",
    "要检查哪些变量是这些变量，我们可以按个体对数据进行分组并获得标准差的总和。如果它为零，则意味着对于任何个人来说，变量都不会随时间变化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c4ed43ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "year            1334.971910\n",
       "black              0.000000\n",
       "exper           1334.971910\n",
       "hisp               0.000000\n",
       "hours         203098.215649\n",
       "married          140.372801\n",
       "educ               0.000000\n",
       "union            106.512445\n",
       "lwage            173.929670\n",
       "expersq        17608.242825\n",
       "occupation       739.222281\n",
       "dtype: float64"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.groupby(\"nr\").std().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d2aa219",
   "metadata": {},
   "source": [
    "对于我们的数据，我们需要删除实体假人，`black`和`hisp`，因为它们对于个人来说是恒定的。 此外，我们需要取消教育。 我们也不会使用职业，因为这可能会调节婚姻对工资的影响（可能是单身男性能够承担更多时间要求更高的职位）。 选择了我们将使用的功能后，是时候估计这个模型了。\n",
    "\n",
    "要运行我们的固定效应模型，首先，让我们获取平均数据。 我们可以通过按个人对所有内容进行分组并取平均值来实现这一点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3d577f59",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>married</th>\n",
       "      <th>expersq</th>\n",
       "      <th>union</th>\n",
       "      <th>hours</th>\n",
       "      <th>lwage</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>nr</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.000</td>\n",
       "      <td>25.5</td>\n",
       "      <td>0.125</td>\n",
       "      <td>2807.625</td>\n",
       "      <td>1.255652</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.000</td>\n",
       "      <td>61.5</td>\n",
       "      <td>0.000</td>\n",
       "      <td>2504.125</td>\n",
       "      <td>1.637786</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>1.000</td>\n",
       "      <td>61.5</td>\n",
       "      <td>0.000</td>\n",
       "      <td>2350.500</td>\n",
       "      <td>2.034387</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>0.125</td>\n",
       "      <td>35.5</td>\n",
       "      <td>0.250</td>\n",
       "      <td>2225.875</td>\n",
       "      <td>1.773664</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>110</th>\n",
       "      <td>0.500</td>\n",
       "      <td>77.5</td>\n",
       "      <td>0.125</td>\n",
       "      <td>2108.000</td>\n",
       "      <td>2.055129</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     married  expersq  union     hours     lwage\n",
       "nr                                              \n",
       "13     0.000     25.5  0.125  2807.625  1.255652\n",
       "17     0.000     61.5  0.000  2504.125  1.637786\n",
       "18     1.000     61.5  0.000  2350.500  2.034387\n",
       "45     0.125     35.5  0.250  2225.875  1.773664\n",
       "110    0.500     77.5  0.125  2108.000  2.055129"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y = \"lwage\"\n",
    "T = \"married\"\n",
    "X = [T, \"expersq\", \"union\", \"hours\"]\n",
    "\n",
    "mean_data = data.groupby(\"nr\")[X+[Y]].mean()\n",
    "mean_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e61ac3e3",
   "metadata": {},
   "source": [
    "为了将数据围绕均值标准化（demean），我们需要将原始数据的索引设置为个体标识符，`nr`。 然后，我们可以简单地从一个数据集中减去对应的数据均值的数据集。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "43e871e0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>married</th>\n",
       "      <th>expersq</th>\n",
       "      <th>union</th>\n",
       "      <th>hours</th>\n",
       "      <th>lwage</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>nr</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-24.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>-135.625</td>\n",
       "      <td>-0.058112</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-21.5</td>\n",
       "      <td>0.875</td>\n",
       "      <td>-487.625</td>\n",
       "      <td>0.597408</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-16.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>132.375</td>\n",
       "      <td>0.088810</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-9.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>152.375</td>\n",
       "      <td>0.177561</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>0.0</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>-0.125</td>\n",
       "      <td>263.375</td>\n",
       "      <td>0.312473</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    married  expersq  union    hours     lwage\n",
       "nr                                            \n",
       "13      0.0    -24.5 -0.125 -135.625 -0.058112\n",
       "13      0.0    -21.5  0.875 -487.625  0.597408\n",
       "13      0.0    -16.5 -0.125  132.375  0.088810\n",
       "13      0.0     -9.5 -0.125  152.375  0.177561\n",
       "13      0.0     -0.5 -0.125  263.375  0.312473"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "demeaned_data = (data\n",
    "               .set_index(\"nr\") # set the index as the person indicator\n",
    "               [X+[Y]]\n",
    "               - mean_data) # subtract the mean data\n",
    "\n",
    "demeaned_data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "432d2b96",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>-9.021e-17</td> <td>    0.005</td> <td>-1.78e-14</td> <td> 1.000</td> <td>   -0.010</td> <td>    0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>   <td>    0.1147</td> <td>    0.017</td> <td>    6.756</td> <td> 0.000</td> <td>    0.081</td> <td>    0.148</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>   <td>    0.0040</td> <td>    0.000</td> <td>   21.958</td> <td> 0.000</td> <td>    0.004</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>     <td>    0.0784</td> <td>    0.018</td> <td>    4.261</td> <td> 0.000</td> <td>    0.042</td> <td>    0.115</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>     <td> -8.46e-05</td> <td> 1.25e-05</td> <td>   -6.744</td> <td> 0.000</td> <td>   -0.000</td> <td>   -6e-05</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = smf.ols(f\"{Y} ~ {'+'.join(X)}\", data=demeaned_data).fit()\n",
    "mod.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "782d0fd4",
   "metadata": {},
   "source": [
    "如果我们相信固定效应消除了所有遗漏的变量偏差，那么这个模型告诉我们婚姻使男人的工资增加了 11%。 这个结果非常显着。 这里的一个细节是，对于固定效应模型，需要对标准误差进行聚类。 因此，我们可以使用库 `linearmodels` 并将参数 `cluster_entity` 设置为 True，而不是手动进行所有估计（这只是出于教学原因）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4cf76fde",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Parameter Estimates</caption>\n",
       "<tr>\n",
       "     <td></td>     <th>Parameter</th> <th>Std. Err.</th> <th>T-stat</th>  <th>P-value</th> <th>Lower CI</th>  <th>Upper CI</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>  <td>0.0040</td>    <td>0.0002</td>   <td>16.552</td>  <td>0.0000</td>   <td>0.0035</td>    <td>0.0044</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>   <td>-8.46e-05</td> <td>2.22e-05</td>  <td>-3.8105</td> <td>0.0001</td>   <td>-0.0001</td> <td>-4.107e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>  <td>0.1147</td>    <td>0.0220</td>   <td>5.2213</td>  <td>0.0000</td>   <td>0.0716</td>    <td>0.1577</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>    <td>0.0784</td>    <td>0.0236</td>   <td>3.3225</td>  <td>0.0009</td>   <td>0.0322</td>    <td>0.1247</td>  \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from linearmodels.panel import PanelOLS\n",
    "mod = PanelOLS.from_formula(\"lwage ~ expersq+union+married+hours+EntityEffects\",\n",
    "                            data=data.set_index([\"nr\", \"year\"]))\n",
    "\n",
    "result = mod.fit(cov_type='clustered', cluster_entity=True)\n",
    "result.summary.tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "621ff0ea",
   "metadata": {},
   "source": [
    "请注意参数估计值与我们使用时间贬值数据得到的参数估计值是如何相同的。 唯一的区别是标准误差有点大。 现在，将其与不考虑数据时间结构的简单 OLS 模型进行比较。 对于这个模型，我们添加了及时不变的变量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a133ca1c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    0.2654</td> <td>    0.065</td> <td>    4.103</td> <td> 0.000</td> <td>    0.139</td> <td>    0.392</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>   <td>    0.0032</td> <td>    0.000</td> <td>   15.750</td> <td> 0.000</td> <td>    0.003</td> <td>    0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>     <td>    0.1829</td> <td>    0.017</td> <td>   10.598</td> <td> 0.000</td> <td>    0.149</td> <td>    0.217</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>   <td>    0.1410</td> <td>    0.016</td> <td>    8.931</td> <td> 0.000</td> <td>    0.110</td> <td>    0.172</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>     <td> -5.32e-05</td> <td> 1.34e-05</td> <td>   -3.978</td> <td> 0.000</td> <td>-7.94e-05</td> <td> -2.7e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>black</th>     <td>   -0.1347</td> <td>    0.024</td> <td>   -5.679</td> <td> 0.000</td> <td>   -0.181</td> <td>   -0.088</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hisp</th>      <td>    0.0132</td> <td>    0.021</td> <td>    0.632</td> <td> 0.528</td> <td>   -0.028</td> <td>    0.054</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>educ</th>      <td>    0.1057</td> <td>    0.005</td> <td>   22.550</td> <td> 0.000</td> <td>    0.097</td> <td>    0.115</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = smf.ols(\"lwage ~ expersq+union+married+hours+black+hisp+educ\", data=data).fit()\n",
    "mod.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bb84cbe",
   "metadata": {},
   "source": [
    "这个模型是说婚姻使男人的工资增加了 14%。 比我们在固定效应模型中发现的效应要大一些。 这表明由于固定的个体因素（如智力和美貌）没有被添加到模型中，结果存在一些省略变量偏差。\n",
    "\n",
    "## 固定效应的可视化\n",
    "\n",
    "为了扩展我们对固定效应模型如何工作的直觉，让我们稍微转向另一个例子。 假设您在一家大型科技公司工作，并且您想估计广告牌营销活动对应用内购买的影响。 当您查看过去的数据时，您会发现营销部门倾向于花费更多的钱在购买水平较低的城市放置广告牌。 这是有道理的吧？ 如果销售额猛增，他们就不需要做很多广告了。 如果您在此数据上运行回归模型，看起来营销成本较高会导致应用内购买减少，但这只是因为营销投资偏向于低支出地区。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2bc0d147",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAE0CAYAAACirQ3aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABe30lEQVR4nO3dd1hUx9fA8e9SbSiINAuoiAXFXhDsBTTYG/YSu8afJXaNiSW2qEmMBns3xti72BWJGhVb1CiKotgFEUSQtu8fvGwktF1c3EXP53nyhL1l7mHc3cOdmTujCA8PVyKEEELkQAa6DkAIIYTIKkliQgghcixJYkIIIXIsSWJCCCFyLEliQgghcixJYkIIIXIsSWJC5/z8/DA3N2fWrFm6DiWF4OBgzM3NGTx4sK5DEWrYuHGjVt5H+vp+FGmTJCa0LjExkXXr1tGiRQtKlChBoUKFcHR0pHbt2gwePJht27bpOkS98/r1a+bPn0+TJk0oXrw41tbWlCtXjm7durFnz550z/Py8sLc3Bw/Pz+1rnP37l1GjBhB1apVsbW1pXDhwri4uNCqVSu+//57Hj58qFY5yV/05ubmVKhQgcTExDSPe/v2Lfb29qpjb9++rVb5QqjLSNcBiE9LYmIiXbp0wdfXl/z589OsWTMKFy7M69evuXfvHrt37+bChQu0b99edU61atX466+/sLS01GHkunP27Fl69OjBixcvKFWqFO3bt6dAgQIEBwdz+PBh9u3bR5MmTVi1ahX58+fP8nVOnTqFt7c30dHRVK9encaNG5M7d24ePnzI9evX+eGHHyhbtizFihVTu0wjIyNCQkI4cuQIHh4eqfbv2LGDiIgIjIyMiI+Pz3LsQqRHkpjQqq1bt+Lr60uFChXYt28fBQoUSLE/JiaGM2fOpNiWJ08eSpcu/THD1Bu3b9+mY8eOvHnzhhkzZjBkyBAMDP5tIAkNDaVv374cOXKEPn36sGXLlhT7NTFixAiio6NZvHgx3bp1SzMWIyPNvhLq16/Pn3/+ydq1a9NMYuvWrcPKyoqSJUty7ty5LMUtREakOVFoVfIXVdeuXVMlMIBcuXLRsGHDFNvS64NIbiq7f/8+y5Ytw9XVFVtbW1xcXJg/fz5KZdKMadu2baNhw4YULlyYUqVKMWbMGGJiYlJd29zcHBcXF16/fs2YMWMoV64cNjY2uLq6snz5clV56oiJieGXX36hfv36FClShMKFC9OgQQNWrVqlUTljx44lMjKSYcOG8dVXX6VKUJaWlmzYsAF7e3uOHj3K9u3b1S77fS9evCAoKIj8+fOnmcAASpcuTcmSJTUqt0CBArRq1QpfX1+ePXuWYt8///zDuXPn6Nq1a4bJ8dSpU3Ts2JESJUpgbW1NpUqVGDduHC9evEjz+KCgIHr16oWDgwOFCxfGw8ODgwcPZhjn69ev+f7776lduzZ2dnYULVqUZs2asXPnTo1+X6F/JIkJrSpYsCCQ1PeiLd988w1z586lWrVq9OjRg9jYWKZPn86sWbNYtGgRI0aMwMnJiV69epE/f36WL1/OxIkT0ywrLi6ONm3acOzYMdq3b0/Pnj0JCwtjzJgxjB8/Xq14IiMj8fLy4ptvvkGpVNK1a1e6detGREQEo0aNYsiQIWqVc//+fU6cOIGpqSkjR45M97h8+fLx1VdfAbB69Wq1yv6vAgUKYGRkRFRUFE+ePMlSGenp1asX8fHx/Pbbbym2r127FoCePXume+7q1atp3bo1/v7+NGvWjKFDh1KkSBGWLl1Kw4YNU/XR3b17lyZNmrBr1y5q1KjBoEGDKFKkCN27d0+37/Dx48c0atSIH374AXNzc3r37k379u0JDg6md+/eMoAjh5PmRKFVLVu25KeffmLVqlVERETQvHlzKleuTIkSJVAoFFkq8/r16/z5559YW1sD8NVXX1GjRg1++eUX8ufPz6lTpyhRogQA48ePp2rVqmzYsIEJEyZgZWWVoqynT59SvHhx/vzzT0xNTVXnNGzYkKVLl9KuXTtq1aqVYTwTJ07k4sWLfPfdd4wYMUK1/d27d/To0YNNmzbRsmVLvvjiiwzLOXv2LACVK1fGwsIiw2MbNWoEwIULF0hISMDQ0DDD4//LxMQELy8vdu3aRbNmzejduzeurq5UqFABMzMzjcr6Lzc3N0qXLs26desYMWIECoWCd+/esXnzZurUqYOjo2Oa5z148IBx48aRJ08ejhw5Qrly5VT7ZsyYwbx58/j666/5448/VNtHjx5NWFgY06dPZ9iwYartvr6+eHt7p3mdwYMHExQUxIoVK+jQoYNqe0REBC1atGDu3Ll4eXlRsWLFD6oHoRtyJya0qmLFiixbtgxra2u2bNnCl19+SdWqVSlevDje3t7s3LlTo+Y2SPriSk5gAA4ODri6uhIdHc2XX36pSmCQdMfRvHlzYmNjuXXrVprlTZkyRZXAIKnJbtSoUUDSMO2MvHr1ik2bNlGxYsUUCQzA1NSUKVOmALB58+ZMf6/k5rciRYpkemzyMe/evSMsLCzT49Py888/07JlSx48eMDUqVNp3rw59vb21K5dmylTphASEpKlcgF69OjBvXv3OHXqFAC7d+8mLCyMXr16pXvOH3/8QWxsLH379k2RwADGjBmDnZ0dhw4d4vHjxwA8evSI48ePU7Ro0VSPPXh6elK3bt1U17h+/TonT57Ey8srRQIDyJ8/P+PHj0epVLJly5Ys/d5C9+ROTGhd27ZtadGiBX5+fpw5c4br169z9uxZfH198fX1xcPDgw0bNmBiYqJWeWn9hWxrawuAi4tLuvuSv/zeZ2RklOadlru7OwBXr17NMJaLFy8SHx+PgYFBms1QySPwAgMDMywHUCVzde5Q3z8mq3e05ubmrF+/nuDgYI4dO8bly5e5cuUKV69e5ebNm6xatYq1a9fSuHFjjcvu2rUr06dPZ926ddSvX5+1a9diYWFBq1at0j3nypUrANSrVy/VPlNTU1xdXdmxYwdXr16lcOHCqn8bV1fXNPvY3N3dUz1qkNxHGxkZmea/V2hoKIAM/c/BJImJbGFsbEyjRo1UzWCJiYns3r2boUOHcujQIVatWsWgQYPUKiut5q7k5rSM9sXFxaXaZ2lpmWZTXHKzY0RERIaxJN8FXb58mcuXL6d73Js3bzIsB8DGxgZArTug5GNMTEwybXrMjIODA3369FG9fvLkCaNGjeLAgQMMGjSIGzduYGxsrFGZlpaWeHl5sXfvXs6fP4+/vz8DBw5Mccf7X8l1/f5d9vuS6yf5uOT//7eJOFla5ST/e508eZKTJ0+mG0tUVFS6+4R+k+ZE8VEYGBjQpk0bVTPQiRMndBJHaGgoCQkJqbYnj4TL7Dms5P0DBgwgPDw83f8yu6MDqF27NpCUEMPDwzM8Nrm+atSooXF/WGbs7OxYuXIlJiYmvHjxguvXr2epnN69e/Pu3Tt69+6NUqnMsCkR/q3L58+fp7k/ubk1+bjk/6c3ajGtcpLPmTFjRob/Xnv37lXjNxT6SJKY+KiS75w07RfTlvj4+DSfV/L39wfSbrp8X/Xq1TEwMEj1rFtWFC9enHr16vHu3Tt++umndI+Liopi0aJFACnuoLQpV65cajfvpqdevXoUL16cR48eUatWrVT9XP9VqVIlgDRnG3n37p3q3yn5uOR/m3PnzqX54HTyv+H7atasCaCVfy+hnySJCa3aunUrx48fT3MaomfPnrFu3Trg3z4oXZg+fTrv3r1TvQ4NDWXBggUA6T5DlaxQoUJ4e3tz7do1Zs2aleaX6aNHj9TuY5kzZw558+Zl4cKF+Pj4pEruYWFh9OjRg/v379O4cWPatWunVrn/FRUVxdy5c9O96/n111958+YN5ubmmSaf9CgUCtatW8eGDRtYuHBhpsd36tQJExMTVq5cmaq+FixYwOPHj/Hw8MDOzg5IGtySPOzex8cnxfG+vr5pJsPKlSvj7u7O/v37Wbt2bZp/PN25c0ft6baE/pE+MaFVFy5cYMmSJaqHiB0cHICkyXQPHTpEdHQ0NWvWpH///jqJz9bWlnfv3uHm5kbz5s159+4du3fv5tmzZwwcODDT4fUAc+fOJSgoiDlz5rB582bc3NywsbHh2bNn3Llzh/Pnz/P999+rNQtJuXLl+OOPP+jZsycTJkxg1apVNGjQADMzMx4+fIivry8RERGqaafSm63jp59+SvWcVrKePXtSrlw5Zs6cyZw5c6hWrRouLi6Ym5vz6tUrzp49y40bNzAyMuLnn3/OsB8rMxUrVlR7qLq9vT1z5sxh1KhRNGzYkDZt2mBjY8O5c+fw9/enSJEizJ8/P8U58+bNo2nTpnzzzTecPHmSihUrcv/+fXbv3k2zZs3SfOh5xYoVtG7dmuHDh7N06VJq1KiBhYUFjx8/5p9//uHq1ats2LBBo+m2hP6QJCa0atiwYTg5OXH8+HFu3LjB8ePHefv2LRYWFtSsWZM2bdrQvXt3jQcOaIuxsTE7duxg+vTpbN26lbCwMEqUKMHXX3+tdmI1MzNj7969rF+/ni1btrB3715iYmKwsrLC3t6eKVOm0KZNG7Vjcnd35+LFiyxfvpyDBw+yefNmoqOjKViwIO7u7nTp0oWWLVtmOCrx6NGj6e6rU6cOtWrVYtu2bRw/fpyzZ89y4MABXrx4gbGxMcWKFaNPnz4MGDAgy3dhWdWnTx9KlizJL7/8wr59+4iKisLOzo4BAwakerQCwNHRkSNHjvDdd99x4sQJ/vzzT8qXL8/GjRt5+fJlmknMzs6O48ePs3z5cnbt2sW2bduIi4vD2tqaUqVKMXv2bOrUqfOxfmWhZYrw8HDddE4I8ZGZm5tTrFgxrl27putQhBBaIn1iQgghcixJYkIIIXIsSWJCCCFyLBnYIT4bmT1QLITIeeROTAghRI4lSUwIIUSOJUlMCCFEjiVJLAdQZ1kP8S+pL81IfWlG6ksz2V1fksSEEELkWJLEhBBC5FiSxIQQQuRYksSEEELkWPKwsxACSFowNCoqStdh6L1cuXLx+vVrXYeRY6hTX3nz5sXIKGvpSJKY+GQER8YxIyCSoJemlHwcxuSqZjiY6WbJl5wmPj6eyMhIzM3NM1zyRYCpqSm5cuXSdRg5Rmb1pVQqCQ8Px8zMLEuJTJKY+CQER8bRxjeUe5EJgCEXI6K58CKWnZ6WksjUEBUVJQlM6IRCocDc3JyIiAgKFCig8flqJ7HQ0FDOnj3L7du3CQ0NRaFQYGlpSenSpalVqxaWlpYaX1wIbZkREPn/Cexf9yITmBEQyfL6BXUUVc4iCUzoyoe89zJMYu/evWPLli1s3LiRc+fOoVSmvX6mQqGgZs2adOvWjU6dOn3Q8uZCZMWTtwlpbn+aznYhxKch3dGJq1evpkqVKowaNYr8+fMzY8YMDhw4wM2bN3n69ClPnjzh5s2b7N+/n2nTpmFmZsbXX39NlSpVWL169cf8HYTALo9hmttt09kuhPg0KMLDw9O8vXJ2dmbw4MH06NEDc3NztQoLDw9n/fr1LFmyhOvXr2szzs9aYGAgTk5Oug5Dr6XsE0tSwsxQ+sTUEBgYiLW1dZb6Iz5HMTExH2VgR3BwMJUqVeL48eNUqVIl26+XnlmzZrF7927OnDmTpfPVra/Xr19n6T2YbhKLi4vD2DhrH/4POVekJklMParRiaFvKGmZT0YnqiknJ7HBgwezadMmAAwNDbGzs8PDw4MpU6ao/ce3pj5WEktISODly5dYWlpmefi5Ovz8/GjZsiV3795Nc2zDmzdviI2NpWDBrPUtZ3cSS7dmPiQJSQITuuBgZszy+gUJDAzFycle1+GIj6RBgwYsXbqU+Ph4bt26xVdffcXr169ZuXJltl43u/9YNzQ0xMbGJtvKV1e+fPl0HUKG1J6xIyoqigsXLrBnzx62bNnCnj17uHDhgjwcKYTQKVNTU2xsbChSpAiNGjWibdu2HDt2LMUxGzZsoFatWtjY2FCtWjUWL15MYmKiav+dO3f44osvsLGxoXr16hw6dIgiRYqwceNGIKlpz9zcnK1bt9K+fXtsbW1Vff+Zlb169WqqVauGjY0Njo6OtGvXjvj4eACuX79Oq1atKFasGEWLFsXd3Z1Tp06luOalS5dUZfn7+9O4cWNsbGxwcnJiwoQJxMbGqvZ7eXnx9ddfM23aNEqWLEmpUqWYPHlying0NWvWLGrXrq16PXjwYLy9vfHx8aFcuXI4ODgwZMgQ3r59qzpGqVTy888/U7lyZYoXL46bmxubN2/OcgwZyfQe9fbt20ybNo3Dhw8TFxeXYoSiQqHA2NiYJk2a8M0331C2bNlsCVIIoRs/O8z7qNcbHjz6g86/f/8+R48eTXGHtHbtWmbOnMncuXOpVKkSN2/eZPjw4RgbGzNgwAASExPp3r071tbWHD58mJiYGCZMmMC7d+9SlT916lSmTJnC4sWLMTY2zrTsS5cuMXr0aHx8fHB1deX169eqJAXQv39/KlSowNGjRzEyMuL69evpNr09fvyYjh074u3tza+//sq9e/f43//+h4GBAd9//73quC1btjBw4EAOHTrEtWvX6NevH5UrV6ZDhw4fVLfvO3PmDDY2NuzcuZNHjx7Ru3dvSpUqxahRowCYMWMGu3btYt68eRQrVoyrV68yfPhwzM3N8fT01FockEkSu3z5Mi1btsTExIRevXpRvXp1bG1tyZUrFzExMTx9+pTz58+zfft2mjZtyp49e6hcubJWAxT6J7nv6cnbBOzyGErfk9CpI0eOUKRIERISEoiJiQFI8aX+ww8/MHXqVFq3bg1A8eLFuXfvHitXrmTAgAEcP36cwMBAtm/fTuHChQGYOXNmml+2AwYMoGXLlqpEk1nZDx8+JG/evDRv3hwzMzMAXFxcVOU9fPiQr776itKlSwNQsmTJdH/PlStXYmNjw/z58zEwMKBMmTJ8++23jBw5kkmTJpEnTx4AypQpw6RJkwAoVaoUa9eu5eTJk1pNYmZmZixYsAAjIyPKlClDmzZtOHnyJKNGjSIqKorFixezfft23NzciImJoUyZMly8eJEVK1Z83CQ2ZcoUihUrxr59+7CwsEjzmE6dOjFhwgS8vLz49ttv2bVrl1YDFPolrVGAMjOG0CU3Nzd+/vlnoqOjWbt2Lffv32fQoEEAvHz5kpCQEEaOHMnXX3+tOic+Pl7VqnT79m3s7OxUCQygatWqGBik7m15f5SgOmU3bNiQokWLUqlSJRo3bkzDhg1p2bKlKqENGTKE//3vf2zatIn69evTqlUrVUL7r1u3blGjRo0UcdWuXZvY2FiCgoKoUKECAOXLl09xnq2tLS9evFCjJtVXpkyZFINNbG1tuXDhgirOmJgYOnTogEKhQKlUolAoiIuLw95e+33VGSaxgIAApk6dmm4CS1awYEG+/PJLpk6dqtXghP6RmTGEvsmTJ4/qDmbu3Lm0aNGCuXPnMmHCBFVf0IIFC6hVq1aa56c3iUNa8ubNq/pZnbLNzMw4deoU/v7+nDhxgh9//JHp06dz7Ngx7OzsmDBhAp06deLw4cMcO3aMOXPmsGDBAnr06JFmnOnNbPH+9v8ONklOJNqU0TWS62XTpk0UK1aMd+/eqSbAyI5RlhmWaGJiQkREhFoFRUZGYmJiopWghP6SmTE+Lx/aR6UL48aNo2PHjvTu3Vt1h3Xv3j26dOmS5vFlypThyZMnPHnyBDs7OwAuXbqU6WAIa2vrTMuGpC/u+vXrU79+fSZMmECpUqXw9fWld+/eADg6OuLo6MigQYMYNWoU69evTzOJlS1blh07dpCYmKi6Gztz5gwmJiaUKFFCnar5KMqUKYOpqSkPHz6kfv362f5IQoZJrEmTJixatAg3N7d0/9IAOHfuHIsWLaJp06ZaD1DoF5kZQ+i7unXrUrZsWebNm8f8+fMZP348Y8eOpUCBAnh4eBAXF8eVK1d48uQJo0aNomHDhjg5OTF48GCmT59OTEwMkyZNwsjIKNM5/TIr++DBg9y7dw83NzcsLCzw8/PjzZs3lC5dmujoaL755htat26Nvb09L1684OzZs1SrVi3Na/Xt2xcfHx++/vprBg0axP3795k6dSr9+/dX9Yd9iBs3bqR6Tiu5iVITZmZmDBs2jG+++QalUkn16tWJi4vjwoULGBgYqJK3tmSYxL7//nuuXr1K8+bNqVChgmqYqKmpKe/evePZs2dcvHiRv//+GycnJ2bMmKHV4IT+mVzVjAsvYlPNjDG5qpkOoxIipaFDhzJ06FCGDx9Oz549yZMnDwsXLmTatGnkypWLcuXK0b9/fwAMDAzYsGEDw4YNo3Hjxtjb2zNjxgx69OiR6R1EZmUXKFCAffv2MXfuXKKjoylRogQLFy7Ezc2N2NhYwsPDGTx4MM+fP6dgwYJ4enoyffr0NK9VuHBhtmzZwpQpU6hbty4FChSgQ4cOTJkyRSt11rJly1TbQkJCslTWpEmTsLKyYtGiRdy7dw8zMzNcXFwYPnz4h4aZSrozdiSLiYlh5cqV7Ny5k7///ls1+geSns9wcXGhdevW9O3bl9y5c2s9QKF/M3Ykj058+jYBWz0cnahv9aXvcvKMHdnl2rVr1K1blxMnTqQacf2xZuz4VOhs2qm0KJVKXr16RXR0NLlz58bCwkKWb/gI5EtZM1JfmpEkBnv27CFv3ryULFmSBw8eMGnSJJRKJX5+fqm+4ySJaUZn006lRaFQZHn+LCGE0Fdv3rzhu+++49GjR5ibm1OnTh1mzpwpf6TnAGolsbi4OE6fPq3qsEzOrHZ2dlSqVIk6derIfIlCiByrS5cuGY4wFPor0yS2efNmpkyZwosXL1TPARgbGxMXFwck3Z0VKlSIqVOnavwm8Pf355dfflElx8WLF9OtW7c0jx0+fDhr165l+vTpDBs2TKPrCCGE+DRlOAHw1q1bGTRoEKVLl2bNmjVcv36dsLAwnj9/TlhYGNevX2f16tWULl2aoUOHsnXrVo0uHhUVhbOzM7Nnz85wUMiuXbsICAhQPcPxuQiOjKP/yTAGXTWl/8kwgiPjdB2SEELolQwHdri7u1O0aFG1Zh/29vYmJCQEf3//LAVSpEgR5s6dm+pO7MGDB3h6erJz5046dOjAgAEDPos7MVnkMetkYIdmkgd25M+fX/qA1CADOzSjTn0plUoiIiKyNLAjwzuxu3fv4uXlpVZBXl5e3L17V+MAMhIfH0+/fv0YPXo0ZcqU0WrZ+i6j6Z2E0La8efMSHh6u9emJhMiMUqkkPDw8xZRemsiwT6xQoULcvHlTrYJu3LhBoUKFshREembNmoWFhQV9+/ZV+5zAwECtxqArQS9NgdSzYASFviEwMPTjB5TDfCrvg4/l3r17AERERKQ58a0Q2SUxMZHY2FhevnyZ5v7MWlUyTGJdunThxx9/pFChQvTt2zfN5b7Dw8NZuXIlK1asUK0low2nT5/mt99+w8/PT6PzPpVmpJKPw7gYEZ16u2U+WbU4E9KcqBmpL81IfWkmu+srwyQ2btw4Hj16xIwZM5g1axYlSpRINe3UvXv3SEhIwNvbm7Fjx2otMD8/P54+fZqiGTEhIYFvv/0WHx8fbty4obVr6SOZ3kkIITKXYRIzMjLi119/ZcCAAezcuZOrV6/y9OlT1YwdRYsWpUWLFrRu3Vrri2H269dPtdBcsvbt29O+fXt69eql1WvpIwczY3Z6WjIjIJKg0DeUtMynd9M7CSGErqn1sHPlypWzZcXmN2/eEBQUBCS1i4aEhHD16lUsLCwoVqwYVlZWKYM1MsLGxuazuZV3MDNmef2CBAaGShOiEEKkQac9uJcuXaJevXrUq1eP6OhoZs2aRb169Zg5c6YuwxJCCJFDqHUndu3aNXbu3JnmtFOVK1emVatWVKxYUeOL161bl/DwcLWPv3btmsbXEEII8enK8E4sPj6eYcOGUb9+fX788UeCgoIoUKAADg4OFChQgKCgIBYsWECDBg0YOnQo8fHxHytuIYQQIuM7sR9++IHffvuNsWPH0r9/fywtLVMdExYWxtKlS5k3bx5FixZlwoQJ2RasEEII8b4M78R+++03+vfvz/jx49NMYAAFCxZkwoQJ9O3bl40bN2ZLkEIIIURaMkxiL168wNnZWa2CKlSokO4T10IIIUR2yDCJOTo6cuDAAbUK2r9/P46OjloJSgghhFBHhkls+PDhHDx4kHbt2nHgwAGePXuWYv+zZ8/Yv38/7dq149ChQwwfPjxbgxVCCCHel+HAjk6dOhEfH8+3336bYomU9xfFVCqVWFpasnDhQjp16pS90QohhBDvyfQ5sa5du9KhQwf8/PxSPSdma2tL5cqVqVOnDqamph8jXiGEEEJFrYedTUxMaNy4MY0bN87ueIQQQgi1ycJBQgghciytJbGDBw8ydOhQbRUnhBBCZEprSezvv/9m06ZN2iruowuOjKP/yTBaHHhB/5NhBEfG6TokIYQQmVCrT+xTFxwZRxvf0BQLUF54EctOT0tZv0sIIfRYhknMzs5O7YISEhIyP0hPzQiITJHAAO5FJjAjIJLl9QvqKCohhBCZyTCJxcXF4eDgQK1atTIt6Pr16zl2qZQnb9NOwE/T2S6EEEI/ZJjEypYtS548efj1118zLWjevHk5NonZ5TFMc7ttOtuFEELohwwHdlSpUoVr16598uuETa5qRgmzlAmrhJkhk6ua6SgiIYQQ6sjwTqxdu3YkJCQQGhqKjY1NhgU1b96cwoULazW4j8XBzJidnpbMCIjk6dsEbPMkJTAZ1CGEEPotwyTWsGFDGjZsqFZB5cuXp3z58loJShcczIxlEIcQQuQwMmOHEEKIHEuSmBBCiBxLkpgQQogcS5KYEEKIHEuSmBBCiBxLkpgQQogcS5KYEEKIHEujWewDAwPZsGED9+/f59WrVyiVyhT7FQoFu3fv1mqAQgghRHrUTmLbtm1j4MCBGBoa4uTkhLm5eapj/pvUhBBCiOykdhKbOXMmzs7ObNu2DSsrq+yMSQghhFCL2n1ijx49omfPnpLAhBBC6A21k1jp0qUJDQ3NzliEEEIIjaidxKZMmcLq1au5c+eO1i7u7+9P586dKVeuHObm5mzcuFG1Ly4ujm+//RY3NzcKFy5MmTJl6NevHw8fPtTa9YUQQuRsaveJHThwACsrK9zc3KhXrx5FixbF0DDlGlwKhYJ58+apffGoqCicnZ3p0qULgwYNSrHv7du3XLlyhdGjR+Pi4kJERASTJ0+mQ4cO+Pv7Y2Sk0cDKTD25+Jh7x4II+fMBDaY1xtol46VnhFBXcGQcMwIiefI2ATtZ5kcIrVKEh4erNaTQwsIi88IUCsLCwrIUSJEiRZg7dy7dunVL95h//vkHV1dX/P39tbrsy8tbL/it+TqUCSmrwrRALjzmN6Nk01Jau1ZWBAYG4uTkpNMYchJ9qq/gyDja+IZyLzJBta2EmSE7PS31JpHpU33lBFJfmsnu+lL7dubVq1fZFoS6IiMjAdIc3v8hnl1+miqBAbx7HcOefjtVr+tNaUilXlUwMJJnxIV6ZgREpkhgAPciE5gRECnr1wmhBdptk8tGsbGxTJ48mWbNmlGkSJF0jwsMDNS47ER7JSbmpsSGv8vwuFPTjnNq2nEAircrSeneZTHK83H+ms7K7/U505f6CnppChim3h76hsBA/RkopS/1lVNIfWnmQ+ors7u4HJHE4uPjGTBgAK9fv2bTpk0ZHpul21YnKHnYkbX1VxD3Nk6tU+5vD+L+9iAASjRxpMG0xuQvkl/za6tBmi80o0/1VfJxGBcjolNvt8yHk5O9DiJKTZ/qKyeQ+tKMzpoTK1asiIGBAefPn8fY2JiKFSuiUCgyLEyhUHD58mWtBhgfH0/fvn25ceMGe/fupWDB7GmCyWudlyE3hwPwJOAxR8f5Enpbvb+U7x25y70jdwGwLFuIJnM9sa1kly1xipxlclUzLryITdUnNrmqmQ6jEuLTkW4Sc3d3R6FQYGBgkOL1xxQXF8eXX37JzZs32bt3LzY2H2fEoF3VwnQ/3AeAiIevOT7lKPePBal1bug/L9ncKulRAeM8xjSd14xSX5T+6HUn9IODmTE7PS2ZERDJ07cJ2MroRCG0Kt0k5uPjk+FrbXjz5g1BQUnJITExkZCQEK5evYqFhQV2dnb06tWLS5cusWnTJhQKBc+ePQMgf/785M6dW+vxpCV/sQK0Xt0OgHeR7zgz7zRX1lxS69y4t3HsH7JH9brOhHpU7lsNQ+PUfSTi0+VgZiyDOITIJmoPsc8Ofn5+tGzZMtX2Ll26MH78eCpVqpTmeYsXL85wKP7HkBifyOXVAfjNOJGl8126V8J9bF1MC+TK9Fhpg9eM1JdmpL40I/WlGZ31id25c4dSpbL2fJS659atW5fw8PB092e0T9cMjAyo2r86VftXB+CubyCHRh8kNiLjEY7Jrm24wrUNVwBwqF+chjOaUMDePLvCFUKIT1K6SczV1ZXWrVvz5Zdf4u7unmlBSqUSPz8/Vq5cyf79+3nx4oVWA9V3jp5ODPZM+mvj2dWnHBnny8sb6tVB8Mn7rKm7AgALx4I0metJ4erpP0YghBAiSbpJ7NixY0yfPp0WLVpga2tL3bp1qVKlCg4ODpibm6NUKgkPDyc4OJjLly9z6tQpnj9/TuPGjTl69OjH/B30jk1FW7od6AVA5OMITn53nLu+6j0n8epuGFvaJz1GYGhqSNMfmqEoJ31oQgiRlkz7xG7evMmGDRvYt28fwcHBSSf9/0i75EUwHRwc+OKLL+jevTvOzs7ZHHLOFfsmlrML/Lm08mKWzq89pg7VBtTA0ESSWkakz0IzUl+akfrSTHbXl0YDO54+fcrt27dV8yMWLFiQMmXKfLSh75+SxIRErq2/zIlvj2Xp/PKdXagzoR65zD/OKM2cRL5kNCP1pRmpL83ozdyJALa2ttja2mZXLJ8VA0MDKvWuSqXeVQEIOnqXI6MPEh2WenaHtFz//RrXf78GQDF3exp+3xSLEplP0iyEEJ+SHDHt1OegZGNHBlwaCsCL6885OuEQz648Vevch/4PWNdgJQAFHMxp+oMnRWoVy7ZYhRBCX0gS00NW5a3pvLs7AG+evWH/mN08OflYrXNfB4eztdPmpBcK8JjfnLLtnGXGECHEJ0mSmJ7LZ5OPKt/UoNM6J+LexnJu4Vku+vyl3slKODTqAIdGHQCg1kg3qg+uiZGp/LMLIT4N8m2WgxjnMaHO+HrUGV8PZaKSa79d4fikI2qff+7HPzn3458AlOtQnrqT6pO7YJ7sClcIIbKdJLEcSmGgoGL3ylTsXhmA+yfucWTMQaKeR6l1/s2t17m59ToARVyL0uj7phQsZZld4QohRLaQJPaJKN6gBP3ODwbg5a0XHBt/mCcB6vWjPTobwvrGqwHIV9iMpj80w76OQ7bFKoQQ2qJREouIiGDZsmWcOnWKly9fsnDhQqpXr05YWBjr16+nRYsWODo6ZlesQk2FyljRaUdXAKJeROE34wS3dt5U69w3jyPZ0W2L6nWTuZ44d6ogA0PEZy84Mo4ZAZEEvTSl5OMwWVJHT6idxB4/fswXX3zBo0ePcHR05Pbt20RFJTVdFSxYkHXr1vH48WPmzJmTbcEKzeW1ykuzn71o9rMX8TFx/PXLOc4vOqv2+UfG+nJkrC8ANb5ypeawWhjlkg+u+LwER8bRxjf0/xc3NeRiRDQXXsSy09NSEpmOqZ3EvvvuOyIiIjh58iQ2NjapZqn38vLi0KFDWg9QaI9RLmPcxtTBbUwdlIlKbmz5W5Wg1HF+0VlVAizTphz1vmlAnkJ5sytcIfTGjIDIFKtzA9yLTGBGQKSsFadjaiexI0eOMHDgQJydnVXTTr2vePHiPH6sXh+M0D2FgYLy3i6U93YB4MHpYA6POcibx5FqnX9r580UTZRf/NoSJ68y2RKrELr25G1CmtufprNdfDxqJ7G3b99mOEfi27dvSUxM1EpQ4uOzr+NA3zMDAQi7E8qxSYd5dDZE7fOTVrBOWsW6av/q1J3cIBuiFEI37PKkPem2bTrbxcejdhJzdHTk4sWL9O7dO839R44ckRnsPxEFS1nSYXNnAN6GvuX0zJOq4fjqCFh+gYDlF4B/Zx8xMDLIlliF+BgmVzXjwovYFE2KJcwMmVzVTIdRCdAgifXq1YtJkybh7u5OkyZNgKQlWaKiopg9ezanTp3Cx8cn2wIVupHHMg8e85vjMb858THxXFjyl+qBaXW8uP6cXxwXqF4PuDREHrAWOY6DmTE7PS2TRieGvqGkZT4ZnagnNFqKZeTIkaxZswYzMzMiIyOxtLQkPDychIQEBg4cyOzZs7Mz1s+WPi79oExUsrPnVh74BWe5jC57e2Dtov1lfPSxvvSZ1JdmpL40o1dLsfz444907tyZHTt2EBQURGJiIiVKlKBdu3bUrl07u2IUekhhoKDtho6q138tPMOZ+f4albGpxXrVz54/fkHZdtIcLYTQjEZ3YkI3ctpffvePB7Gr9/Ysn+/SvRKNvm+a5fNzWn3pmtSXZqS+NKM3d2LR0dFERUVRqFAh1baXL1+ybt06wsPDad26NdWqVcuWIEXOUrxhSYYHjwYg/P4r1tZfqdH51zZc4dqGKwBYOBakm28vDI1lFJgQIjW1k9jIkSO5efMmJ0+eBCAqKopGjRrx8OFDAHx8fNizZw+urq7ZE6nIkcyLW6gSWmxULCtq+BAXFaf2+a/uhrGo1I+q1/0uDCavlTxgLYRIonYSO3v2LJ07d1a93rp1Kw8fPmTr1q24uLjQrl075s2bx9atW7MlUJHzmeQ1YciN4QAolUr2DthF0KE7GpWxovq/I2C9d3bDtoqdVmMUQuQsaiexZ8+eUaRIEdXrAwcOULNmTRo3bgxAt27d+PHHH9M7XYgUFAoFLZe3Ub2+uPQ8p2ee1KiMzW02qn5uPNuDCl0qais8IUQOoXYSy5s3L+Hh4QDEx8fz559/MnjwYNX+3LlzExmp3pRFQvxXtYE1qDawBpA0Bdb7M+mr4+j4QxwdnzR3ZxHPYjgtk453IT4HaiexKlWqsH79eurVq8eBAwd48+YNzZo1U+2/d+8e1tbW2RKk+LzY13FQ9aNFhLxmtftyjc5/5PuQnx3mqV4PvTUCo1yydJ4QnyK1P9mTJ0+mbdu2NGzYEKVSSatWrahSpYpq/969e6lVq1a2BCk+X/mLFlAltPiYOFa5Lyf65VuNylhc5ifVz739+lHA3lyLEQohdEntJFapUiXOnz/PuXPnMDMzo27duqp94eHh9OvXD3d392wJUghIWkpmwMUhQNLAEN/h+7m1S73FPpOtqbtC9XOrVW0p0VgWcRUiJ5OHnXMAebgyc1fWXuLElKNZPr/6kJq4j6unxYhyDnl/aUbqSzN687Dz+yIjI4mIiEhz6ZVixYp9cFBCaKpSrypU6pXUvO2/4TQXJqm/ejXAhV//4sKvfwFgWaYQ3Xx7oVAotB6nEEK7NEpi69atY+HChQQFBaV7TFoLZqbH39+fX375hStXrvDkyRMWL15Mt27dVPuVSiWzZ89m7dq1hIeHU61aNebNm0e5cuU0CVt8Zqxr2aj60V4/CE/RhKiO0FsvWVh8vur1kJv/wziPiVZjFEJoh9qLPK1fv57hw4dTrFgxJk+ejFKpZPDgwYwcORJra2tcXFz45ZdfNLp4VFQUzs7OzJ49m9y5c6fa//PPP7N48WLmzJnDsWPHsLKyom3btjKUX6itgL05w4NHMzx4NENvjchSGb+WW8jPDvP42WEer+6q/0eaECL7qX0n5uPjQ926ddmxYwdhYWFMnz4dDw8P6tevz7Bhw6hfvz4REREaXdzDwwMPDw8AhgwZkmKfUqnEx8eHESNG0Lp1a1UMTk5ObN26lT59+mh0LSGMchmp7tAA/mj3G08uPtaojHWNVql+/uLXljh5ldFafEIIzal9JxYUFESLFi2STjJIOi0uLmkOPHNzc3r27MmKFZo122QkODiYZ8+e0ahRI9W23Llz4+bmxrlz57R2HfH56rS9q+ourfboOhqfv3/IHtUd2snvjmVDhEKIzGg0Y4dSmTSQMV++fBgaGvL06VPV/oIFC/L4sWZ/1Wbk2bNnAFhZWaXYbmVlxZMnT7R2HSEAag5zpeawpMmrszJjyOXVAVxeHQCAWdH89DndXwaGCPERqJ3EnJycuHHjRtJJRka4uLjw+++/4+3tTUJCAps3b8bBwUHrAf73i0CpVGb45RAYGKj1GPTBp/p7ZZcPqi8b+OJIUhN29Itojnc5pNHpkSERKQaGNN31BcZ59XsZe3l/aUbqSzMfUl+ZDc9XO4l5eXnh4+NDTEwMuXLlYvTo0fTo0YPixYujUCiIiopiyZIlWQ70v2xskpatf/78OUWLFlVtf/nyZaq7s/d9is9vyHMpmtFqfTlBxeCkiYUTYhNY5KT5JNeHW+9X/dzrZF/Mi1toJzYtkfeXZqS+NKM3z4kNGzaMYcOGqV57eXmxf/9+du3ahaGhIc2aNaNOHc37FdLj4OCAjY0Nx48fp2rVqgDExMRw5swZpk2bprXrCKEuQxPDFANDdnTfwgO/YI3KeH+B0NZr2lG8YUmtxSfE5+iDZkV1dXX9oEUw37x5o3rmLDExkZCQEK5evYqFhQXFihVj8ODBzJ8/HycnJ0qVKsW8efPImzcvHTp0+JCwhdCKths6qn4OWHYev+81W0pmV+/tqp9rf+1Ozf/V1lpsQnwudDrtlJ+fHy1btky1vUuXLvj4+Kgedl6zZk2Kh52dnZ11EK3uSPOFZnRdX4/+CmFrx9+zfH4xd3vabuiIwuDjDAzRdX3lNFJfmsnu+tIoia1bt461a9dy//59Xr16lbowhYLQ0FCtBijkQ6MpfaqvqOdRrKjhk/mBGRj09zBMzUy1FFFq+lRfOYHUl2b0pk9s2rRp/PTTT5QvX56OHTtibm6ebUEJ8anIa51X1Y+WmJDIH21/49mVp5mcldKSCv/OhNPjaB8KlrLUaoxC5GRqJ7ENGzbwxRdfsGHDhuyMR4hPloGhAZ13d1e99vv+BAHLLmhUxvrGq1U/t1jWGkfPT++OIDgyjhkBkTx5m4BdHkMmVzXDwUy/H1EQuqN2EouKiqJJkybZGYsQn5W6kxpQd1IDAAL33WL/kD0anb93wC7VzzWG1sJtbN0Mjs4ZgiPjaOMbyr3IBNW2Cy9i2elpKYlMpEntJObq6sr169ezMxYhPltOXmUYHpw0D+PLWy/Y6LFWo/PPLz7H+cVJ07HZVrGj47YuGBiqPauc3pgREJkigQHci0xgRkAky+sX1FFUQp+pncR++OEH2rRpQ6VKlejWrZtMqSNENilUxkrVj/budQxLKi7S6Pynl57wS8kFqtcDr35FrgK5tBpjdnnyNiHN7U/T2S5EukmsVq1aqbbFxsbyv//9j7Fjx1K4cGEMDQ1T7FcoFJw9q9lihEKI9JkWyKVKaMpEJds6b+bRuRCNylj6XhLs5tuLQmXTn/FG1+zyGKa53Tad7UKkm8QKFSqU6m7LysqKUqVKZXtQQojUFAYKOvzRWfX6zLzT/PWLZn80bvT8t5my+aIWlG5ZVmvxacPkqmZceBGbokmxhFnS4A4h0qLTh52FeuS5FM18jvUVdPgOe/rtzPL5VfpWo96UhtoL6AMkj058+jYBWz0cnfg5vr8+hF497Cx0Qz40mvnc6+vV3bAUi3dqqpCzFV329MDAKOcNDPkYPvf3l6b05mHndevWcfjwYdavX5/m/p49e9KsWTO6du2qteCEEJqzcCyo6keLfRPLsiqLSYhVf2DEyxsv+MXx34Eh/QOGkMcyj9bjFEIb1E5iq1atonr16unut7W1ZcWKFZLEhNAjJvlM+CpwJJC0Ft+uXtsIPnlfozKWV/1V9XPnvd2xcbHVZohCfBC1k9jdu3fp1atXuvvLlSvH779nfdJTIUT2UigUtFn37woQfy06y5kfTmtUxu8t/p2xx2NBc8q1L6+1+ITICrWTWGaT+4aFhZGYmKiVoITIiuQBAUEvTSn5OEzvBgTom5pfuVLzK1cCAwMxeWzEzp7bNDr/0KgDHBp1AIAK3SrReGbT7AhT5FAf6/Oo9sCOli1b8vz5c06ePEmuXCkfnIyOjqZhw4YULFiQ/fv3p1OCyCrpSM5cWtMVlTAzlOmK1PDf99frB+Gsqbsiy+UVKG5Oj8N9MDT5NJ/tks9j5j7m51HtO7FRo0bRvn17PD09GTVqFM7OzigUCq5fv86PP/5IYGAgmzdv1mpwQqhLpivSngL25qqBIXFvY1lRaymxEe/UPv/1/XAWOf2oet3v/GDyWufVepxCf33Mz6PaSaxhw4b8+uuvjB07lj59+qi2K5VKzMzM+OWXX2SCYKEzMl1R9jDOY8Lga8OApM/6voG7uesbqFEZ76+n1mlHV+yqFtZqjEL/fMzPo9pJDKBz5854eXlx7Ngx7t+/j1KppESJEjRq1AgzM3miXuiOTFeU/RQKBS2WtVa9DlhxAb/pJzQq44+2v6l+bjSrKS5dK2kpOqFPPubnUa0+sejoaGrXrs2gQYMYNGiQ1oMQGZM2+MxJn1jWaeP99fDPB2zv8keWzy/bzhmPBc1zxMTi8nnMnN71ieXOnZuIiAhMTEy0enEhtMXBzJidnpZJo6FC31DSMp+MTvyIirnZq/rRIh5FsNptmUbn/7P9Bv9svwFAPtt89DrZD6NcGjUUCT3yMT+Pao9OHDRoEOHh4fIsmA7IX36akfrSTHbWV3xMPGvqLifqeVSWy/jy7EDM7PSnu0LeX5rRm2mnRo4cSZ8+fejduzd9+vShRIkS5M6dO9VxVlb6u8yDEOLjMsplRL/zg4GkgSG+I/Zza+dNjcpY5bpU9XOHP7wpUquYVmMUOZvad2IWFhb/npRBu3VYWNiHRyVSkL/8NCP1pRld1dfV9Zc5PvlIls+v/10jKvepqsWI1CPvL83ozZ3Y2LFjc0SnqxAiZ6jYozIVe1QG4PH5R2zpsEmj809+d4yT3x0DoFRzJ77waSXfUZ8htZPYhAkTsjMOIcRnrHCNIqqBIW+evWFlzSUanX/nQCALi88HklbD7ntuIMa5ZVDP50CG/wgh9Eo+m3yqhJYQm8D6Jqt5HRyu9vnvXsfwa9mfVa/7nO5P/mIFtB2m0BNqJ7E5c+ZkeoxCoWDs2LEfFJAQQiQzNDGk96l+qtdHxvpyffM1jcpYXWe56uc2GzrgULe4tsITeiBLAztSFaJQoFQqUSgUMrAjG0hHsmakvjSTU+vr+uZrHBnrm+Xz3SfUo/qgmhqfl1PrS1f0ZmDHq1evUm1LTEzkwYMHLF26lHPnzrF161atBieEEOkp7+1CeW8XAJ5eecLmVhs1Ot9/1in8Z50CwKFBCVqvaScDQ3Igte/EMtOnTx+MjY1ZtkyzJ/VF5uQvP81IfWnmU6uvty+jWF7NJ/MD02FoasiAS0MxyZv2DEWfWn1lN725E8tM3bp1mTp1qraKE0KILMlTKO+/A0PiEtjktZ7QWy/VPj/hXQI+zgtVr3ue6ItFifS7U4RuaS2JBQYGolRq5aZOCCG0wtDYkO6Heqten/j2KFfWXNKojHUNVqp+brWqLdhrKzqhDWonMX9//zS3v379Gj8/P5YvX06bNm20FZcQQmhdg6mNaTC1MQD/7LiB7wjNVqLf/eUO1c/tf+9E0dqS0XRN7STWokWLNDs9lUolhoaGtG/fXq1h+JpISEhg1qxZ/PHHHzx79gwbGxs6derE+PHjMTKSR9yEEFlXtq0zZds6A/Di+nN++2KdRudv6/zv0jNN5zWjXIfyMjBEB9TOBLt37071D6RQKDA3N8fe3j5bFsX86aefWLFiBT4+Pjg7O3P9+nUGDx6MiYmJPI8mxAcKjoxLWirjpSklH4d91kvXWJW3VvWjRb+KZlnlxRqdf3j0QQ6PPghAzeG1qTGkliwl85FobXRidvD29sbCwoIlS/6dgmbQoEG8evWKzZs36zCyj0tGQ2lG6itzsoioehITEtnSfhNPLz3J0vnl2jtTZ1ID8ljm0XJkOUd2fx4NMjtg3bp11KpVCxsbG8qVK8fEiROJjY3NtoDe5+rqyunTp7l9+zYA//zzD35+fjRt2vSjXF+IT9WMgMgUCQzgXmQCMwIidRSRfjIwNMB7ZzeGB49mePBo6kysr9H5N7fdYHnVX/nZYR5bOmwiLDA0myL9fGV4J7ZlyxYGDBhA3rx5KV26NI8ePeLFixcMGDCA2bNnZ3twSqWSGTNmsGDBAgwNDYmPj2f06NFMnjw53XMCAwOzPS4hcrpBV025GGGYanu1AgkscXmng4hynsjgCP7+6SqvrmmemHJZ5aLi6CoUqmadDZF9WjK7i8swiTVt2pSwsDAOHDiAtbU18fHxDBw4kH379hEUFESePNl7i7xt2zamTJnCtGnTKFu2LNeuXWP8+PFMmzaNnj17Zuu19Yk0j2lG6itz/U+GsSUoOtX2jiVzs7x+QR1ElHOk9f56G/qW09+f4Oa2G1kqs/FsD8p7u6Aw+PQGhmT35zHDJFa0aFHGjRvHsGHDVNv+/vtv6tati7+/P87OztkWGED58uX56quvGDx4sGrbDz/8wG+//calS5o965GTyZeyZqS+Mid9YlmX2fsrPiaeCz7nOPfTmSyVX2NoLWr+zxWjXJ/Gv4NOZ+yIiorC1tY2xbbChQsD8PKl+k/AZ9Xbt28xNEzZ5GFoaEhiYmK2X1uIT5mDmTE7PS2TRieGvqGkZb7PenSiNhnlMsJ1pDuuI91RKpXc3HpdNXJRHecXn+P84nMAlG5VlnrfNCSvdd7sCjfHy3QMqC6fe2jWrBk//fQTDg4OlC1blqtXr7J48WI6d+6ss5iE+FQ4mBmzvH5BAgNDcXKSh3azg0KhwLljBZw7VgAg5MwDDo/xJeLha7XOv737H27v/gcA2yp2NJ7tQaGyVtkWb06UYXOihYUFFSpUSHE3FhcXx8mTJ6latSoFC6ZsO1coFPzxxx//LSbLIiMj+f7779m7dy8vX77ExsaG9u3bM3bsWHLlyqW16+g7aR7TjNSXZqS+NKOt+noVFMaxSUcI+fOBxufmLpQHj3nNKN6w5AfHkd102ifm4uKi0Z2YQqHgypUrWglM/Eu+ZDQj9aUZqS/NZEd9Rb+Kxn/WKY0X/EzWcEYTXLpV0suBITpNYkI/yJeMZqS+NCP1pZnsrq/4d/EELD3Pmflpz1ebmaoDquM60g3jPGkvJfOx5ZilWITIiuSpj568TcAuj6EMLhCfPSNTI2r+rzY1/1cbpVLJrV3/cHj0ARLj1BvQFrDsAgHLLgBQ6ovS1P+uEfls8mVnyDolSUzoTFrDvC+8iJVh3kL8P4VCQdk25SjbphwAj86HcGSsL+FBr9Q6/87+29zZnzTjkXUFGxrP8cC6gk22xasLksSEzmQ09ZE8cCtEakVqFKXX8b4AhAeHc+KbIwSfvK/Wuc//fsYmr/UA5DLPRdP5zSnZxDG7Qv1oJIkJnXnyNiHN7U/T2S6E+Je5gzlt1nUAIOZ1DP5z/Ph7o3oD62LCY9jT99+10ep/14iKPStjYJjpdLp6R5KY0Bm7PKnn7gOwTWe7ECJtuQrkovHMpjSe2ZSEuAQurbyI/6xTap9/8rtjnPzuGACV+1aj9ih3TPLpx8CQzEgSEzozuaoZF17Eppr6aHJV7a9NJ8TnwtDYkOqDalJ9UE2USiWB+25x+OuDxMfEq3X+5ZUXubzyIgAlPUrRYGojzArnz86QP0iWktjr168JCQkBkuZXLFCggFaDEp+H96c+evo2AVsZnSiEVikUCkq3KEvpFmUBeBLwmCNjfdVeEibo0B2CDt0BoFA5K5rM9cSmom0mZ31cGiWxM2fOMHXqVP76668U22vWrMmUKVNwc3PTanDi05c89ZEQIvvZVS1MjyN9AIh4+JrjU45y/1iQWue+vPmC31tuAMAknwlN5zXDsZmTTqcmBA2S2KFDh+jWrRv58uWjb9++lCpVCqVSyd27d9m6dSutW7dm48aNeHh4ZGe8QgghtCB/sQK0Xt0OgHcR7zgz7zRX1qq3Okjsm1j2Ddqtel13Un0qf1kNA6OPPzBE7Rk73N3diYuLw9fXFwsLixT7wsLC8PDwwNTUFH//rD1lLtInMypoRupLM1JfmvnU6ysxPpHLqwPwm3EiS+dX7FkZtzF1Mc1vCujRjB13797lm2++SZXAAAoWLEivXr34/vvvtRqc0E8yy4YQnybVZ7uoA3brvmRSlXzE/xnM4dEHiI2MVauMq+suc3XdZQCKNypJ8S8dIRtzvtpJrHjx4kRFRaW7PyoqCgcHB60EJfSXzLIhxKcp/c92cQb//T8Anl19ypFxvry88UKtMu8fC+L+sSCU3yVQuU/VbIlb7QbMcePGsWTJEi5cuJBq3/nz51m+fDkTJkzQanBC/2Q0y4YQIudS57NtU9GWbgd6MTx4NF+eGYCjZym1yj41/TjRYW+1Gm8yte/ETp8+ja2tLR4eHlSpUgVHx6TpSu7evculS5coV64cfn5++Pn5qc5RKBTMmzdP+1ELnZFZNoT4NGn62TYrnJ8Wy9oASQM9zi7w59L/P1/2XwoDBQZG2TOJgdpJbNWqVaqfAwICCAgISLH/xo0b3LhxI8U2SWKfHpllQ4hP04d8tk3ymVBvSkPqTWlIYkIi19Zf5sS3STOA5CmSlzqj66kGemib2kns1Sv1Zk0WnzaZZUOIT5O2PtsGhgZU6l2VSr2T+sD0ZnSiECCzbAjxqcqpn22Nk9jdu3c5fPgwDx8+BKBYsWI0adKEUqXU6+ATOZ/MsiHEpyknfrbVTmLx8fGMGTOGdevWkZiYcoXRiRMn0qNHD+bPn4+RkdzcCSGE+DjUHmI/ffp01qxZQ6dOnTh+/DgPHjzgwYMHHDt2jE6dOrFu3TqmT5+enbEKIYQQKah927Rp0ybatGmDj49Piu1VqlRhyZIlREdH89tvvzF16lStBymEEEKkRe07sbdv31KnTp1099erV4+YmBitBCWEEEI9wZFx9D8ZRosDL+h/MozgyDhdh/RRqZ3E3NzcOHv2bLr7z549K0uxCCHER5Q8VdSWoGhOP41lS1A0bXxDP6tEpnYSmz9/PlevXuXrr7/m1q1bxMXFERcXx61btxg1ahTXrl1j/vz52RmrEEKI98g0cBr0idWoUQOlUsnt27dZvXq1aiE0pTJpJRcjIyNq1KiR4hyFQsHjx4+1GK4QQohkMg2cBkmsbdu2Ol/BUwghxL9kGjgNkth/RyUKIYTQLZkGTqadEkKIHCunThWlTRonsSdPnnDlyhVev36dauYOgC5dumglMCGEEJnLiVNFaZPaSSw2NpavvvqKbdu2kZiYiEKhUA3qeL+vTJKYEEKIj0XtIfYzZ85k27ZtTJgwgb1796JUKvHx8WHHjh00atQIFxcX/P39tR7g06dPGTRoEI6OjtjY2FCrVi1Onz6t9esIIYTIedROYtu2bcPb25vRo0dTrlw5AOzs7GjQoAFbtmwhT548KRbO1Ibw8HA8PT1RKpX88ccfnDt3jrlz52JlZaXV6wghhMiZ1G5OfP78ObVq1Uo66f9nqk+eZkqhUNC6dWt+/PFHfvjhB60Ft3DhQmxtbVm6dKlqW/HixbVWvhBCiJxN7TsxS0tLwsPDATAzMyN37tzcv39ftT8uLo6oqCitBrdv3z6qVatGnz59KFWqFHXq1GHZsmWqvjghhBCfN0V4eLhaGcHb2xtjY2M2bNgAQMeOHbl9+zZLliwhMTGRQYMGUbRoUQ4cOKC14GxsbAAYMmQIbdq04dq1a4wbN45vv/2WAQMGpHlOYGCg1q4vhBBCt5ycnDLcr3YSO3DgABs2bGDlypXkypWLf/75hxYtWhAWFoZSqaRgwYJs2bKFqlWraiVwACsrK6pUqcKhQ4dU26ZNm8bevXv566+/tHYdfRcYGJjpP6T4l9SXZqS+NCP1pZnsri+1+8SaN29O8+bNVa/Lli1LQEAAfn5+GBoa4urqirm5uVaDs7GxoUyZMim2lS5dmpCQEK1eRwghRM70QTN25M+fHy8vL23Fkoqrqyt37txJse3OnTsUK1Ys264phBAi51B7YMd/hYeH07JlS65cuaLNeFIYMmQI58+fZ968eQQFBbFz506WLVtGv379su2aQgghco4sJ7HY2FhOnz6tGrGYHapWrcrGjRvZsWMHtWvXZvr06UycOFGSmBBCCCAHTADs6emJp6enrsMQQgihh7J8JyaEEELoWpaTWK5cuejSpQt2dnbajEcIIYRQW5abE/Pnz8+vv/6qzViEEEIIjWQpiUVFRfHq1as0p3+S4e9CCCE+Fo3WE5s7dy5r164lNDQ03ePCwsK0EpgQQgiRGbWT2NixY1m3bh3NmzfH3d1d67NzCCGEEJpSO4nt3LmTrl27smjRouyMRwghhFCb2qMTExMTqV69enbGIoQQQmhE7STWuHFjzp49m52xCCGEEBpRO4nNnTuXv//+m5kzZ/Ls2bPsjEkIIYRQi9p9Yi4uLiiVSubNm8e8efMwNjbGwCBlDlQoFDx+/FjrQQohhBBpUTuJtW3bFoVCkZ2xCCGEEBpRO4n5+PhkZxxCfLDgyDhmBEQS9NKUko/DmFzVDAczY12HJT4R8v7ST3o/i70Q6giOjKONbyj3IhMAQy5GRHPhRSw7PS3li0Z8MHl/6a8Mk9jFixc1LrBatWpZDkaIrJoREPn/XzD/uheZwIyASJbXL6ijqMSnQt5f+ivDJNakSRO1+8GUSiUKhUKmnRI68eRtQprbn6azXQhNyPtLf2WYxBYvXvyx4hDig9jlMUxzu20624XQhLy/9FeGSaxr164fKw4hPsjkqmZceBGbosmnhJkhk6ua6TAq8amQ95f+kpWdxSfBwcyYnZ6WdCyZm2oFEuhYMrd0ugutkfeX/pLRieKT4WBmzPL6BQkMDMXJyV7X4YhPjLy/9JPciQkhhMixJIkJIYTIsSSJCSGEyLEkiQkhhMixFOHh4UpdByGEEEJkhdyJCSGEyLEkiQkhhMixJIkJIYTIsSSJCSGEyLEkiQkhhMixJInpqQULFtCwYUOKFSuGo6Mj3t7e3LhxQ9dh5Rjz58/H3NycMWPG6DoUvfX06VMGDRqEo6MjNjY21KpVi9OnT+s6LL2UkJDAjBkzqFixIjY2NlSsWJEZM2YQHx+v69D0gr+/P507d6ZcuXKYm5uzcePGFPuVSiWzZs2ibNmy2Nra4uXlxc2bN7VybUlieur06dP07dsXX19fdu/ejZGREW3atOHVq1e6Dk3vnT9/nrVr11K+fHldh6K3wsPD8fT0RKlU8scff3Du3Dnmzp2LlZWVrkPTSz/99BMrVqxgzpw5/PXXX8yePZvly5ezYMECXYemF6KionB2dmb27Nnkzp071f6ff/6ZxYsXM2fOHI4dO4aVlRVt27YlMjLyg68tz4nlEG/evMHe3p6NGzfSvHlzXYejt16/fk39+vX5+eefmTt3Ls7Ozvzwww+6DkvvTJs2DX9/f3x9fXUdSo7g7e2NhYUFS5YsUW0bNGgQr169YvPmzTqMTP8UKVKEuXPn0q1bNyDpLqxs2bL079+f0aNHAxAdHY2TkxPTp0+nT58+H3Q9uRPLId68eUNiYiLm5ua6DkWvjRgxgtatW1O/fn1dh6LX9u3bR7Vq1ejTpw+lSpWiTp06LFu2DKVS/qZNi6urK6dPn+b27dsA/PPPP/j5+dG0aVMdR6b/goODefbsGY0aNVJty507N25ubpw7d+6Dy5elWHKI8ePH4+LiQs2aNXUdit5au3YtQUFBLF26VNeh6L379++zcuVKhgwZwogRI7h27Rrjxo0DYMCAATqOTv+MGDGCN2/eUKtWLQwNDYmPj2f06NH069dP16HpvWfPngGkaqq2srLiyZMnH1y+JLEcYOLEiZw9e5aDBw9iaCjLoaclMDCQadOmceDAAUxMTHQdjt5LTEykSpUqfPvttwBUqlSJoKAgVqxYIUksDdu3b+f3339nxYoVlC1blmvXrjF+/Hjs7e3p2bOnrsPLERQKRYrXSqUy1baskCSm5yZMmMD27dvZs2cPxYsX13U4euuvv/4iNDSU2rVrq7YlJCTw559/smrVKh4/foypqakOI9QvNjY2lClTJsW20qVLExISoqOI9NuUKVP46quvaN++PQDly5fn4cOH/Pjjj5LEMmFjYwPA8+fPKVq0qGr7y5cvtTKQSJKYHhs3bhzbt29n7969lC5dWtfh6DUvLy+qVKmSYtvQoUNxdHRk1KhRcnf2H66urty5cyfFtjt37lCsWDEdRaTf3r59m6oVxNDQkMTERB1FlHM4ODhgY2PD8ePHqVq1KgAxMTGcOXOGadOmfXD5ksT01OjRo9m8eTMbNmzA3Nxc1a6cN29e8uXLp+Po9I+5uXmqQS958uTBwsICZ2dn3QSlx4YMGYKHhwfz5s2jXbt2XL16lWXLlvHNN9/oOjS91KxZM3766SccHBwoW7YsV69eZfHixXTu3FnXoemFN2/eEBQUBCQ1VYeEhHD16lUsLCwoVqwYgwcPZv78+Tg5OVGqVCnmzZtH3rx56dChwwdfW4bY66n0RiGOGzeOCRMmfNxgcigvLy8ZYp8BX19fpk2bxp07dyhatCj9+/dn4MCBWumn+NRERkby/fffs3fvXl6+fImNjQ3t27dn7Nix5MqVS9fh6Zyfnx8tW7ZMtb1Lly74+PigVCqZPXs2a9asITw8nGrVqjFv3jyt/IEpSUwIIUSOJc+JCSGEyLEkiQkhhMixJIkJIYTIsSSJCSGEyLEkiQkhhMixJIkJIYTIsSSJCb3h5eVFjRo1sv06s2bNktUAtGz8+PF4eHik2GZubs6sWbN0FJFuHTx4kKJFi/Ly5Utdh/LJkyQmVDZu3Kia+eLUqVNpHtOoUSPMzc0/SrL5EI8ePWLWrFlcvXpV16Gk6d27dyxduhQPDw/s7e2xsrKiYsWKDB06lMuXL2fLNc+cOcOsWbMIDw/XarkhISGsXr1atVZUdpo/fz5du3albNmymJubM3LkyHSPjYyMZMyYMZQuXRpbW1uaNGnCsWPH0jw2KCiILl26YG9vT7FixejSpQv37t1L89hjx47RpEkTbG1tKV26NGPHjuXNmzcpjvH09MTe3l4WzfwIJImJVHLlysWWLVtSbb979y4BAQE5YoaCx48fM2fOHK5du5Zq35gxY3j69KkOokoSGhpK8+bNGTduHPnz52fcuHHMnz+fjh07cvbsWRo2bMijR4+0ft2zZ88yZ84cXr9+rdVyly5dirW1daq1tZ4+fcqYMWO0eq3p06dz/vx5KlWqlOFxSqWSbt26sX79enr06MGsWbMwMDCgY8eOnDx5MsWxz549o3nz5ly7do2xY8cybtw4rl27RvPmzXn+/HmKY0+ePEnHjh0xNDRk1qxZ9OjRg3Xr1tG1a9cUa7EpFAp69+7NmjVriIiI0F4FiFRk7kSRioeHB7t27WLevHkpZn7fvHkz1tbWODo6arWZ5O3bt+TJk0dr5WXGyMgIIyPdvfUHDx7MpUuXWLVqFe3atUuxb+LEiSxatEhHkWkuPj6ezZs306VLl1TTVWXHHzuXL19WreaQUZPw3r17OXXqFEuWLFHNb9i1a1fc3d2ZNGkSp0+fVh27YMECXr16xblz5yhRogQAzZs3x9XVlQULFjB79mzVsRMnTqRUqVLs2bNHNal0yZIlGTp0KPv27aNFixaqY1u3bs348ePZuXOnzHSfjeROTKTSvn173rx5w8GDB1Ns37p1K+3atcPAIPXbZuPGjbRu3ZrSpUtjbW1NtWrV+Omnn1LN8p3c7/X333/TsmVLChcuzNdff51uLGfOnMHe3p6OHTsSExMDJN0Rfvnllzg6OmJtbY2bmxsbNmxQnfP+irtDhw5VNZEm98+k1Sfm4uJC+/btuXjxIs2aNcPW1pby5cvz66+/poopJCSE7t27U6RIEUqUKMGwYcP4+++/MTc3Z+PGjRnULFy4cIFDhw7RvXv3VAkMkmZGHz58OEWKFFFtu3HjBp07d8be3h47OzuaNm3K4cOHU527YsUK3NzcKFy4MMWLF6d+/fqsWrVK9TtPnToVSFo7LLlO/Pz8gKTk0LFjRxwdHbG1taVSpUoMHDiQqKioDH+fc+fO8fz5cxo0aJBq33/7xJKbq//880+mTZtGmTJlsLW1pW3btty/fz/D6yRTdzmiHTt2YGFhQceOHVXbTE1N6dWrF3///TeBgYGq7Tt37qRJkyaqBAbg6OhIo0aN2LFjh2rb7du3uX79Or169UqxKoK3tzcFChRIcSwkLUFSrlw59u7dq1bMImvkTkykUrhwYdzd3dmyZQutW7cGkr58g4KC6NSpU5pNdMuXL8fJyYkmTZqQO3dujh8/znfffUdERARTpkxJcezr169p164dLVu2pH379hQoUCDNOE6ePEmXLl1o1KgRq1atwsTEhFu3buHp6YmlpSVDhw6lQIECHDp0iK+++oqIiAiGDBlCmTJlGD9+PLNnz6Z3796qNcbKly+f4e8dHBxM586d6dq1Kx07dmT79u1MnDiRsmXLqpZWf/v2La1atSIkJIQBAwZgb2/P3r17GTx4sFp1e+DAAQC1Zz+/c+cOzZo1w8TEhCFDhpA3b15+++03vL29Wbt2rWrS1XXr1jF69GhatWpF//79iYuL459//uHs2bN8+eWXtGzZksDAQLZv387MmTOxtLQEoEyZMrx8+ZK2bdtiaWnJ8OHDMTc3JyQkhAMHDhAVFUXevHnTje/s2bMAVK5cWa3fB5LuZnLnzs3IkSMJDQ1l0aJFDBgwgEOHDqldRmauXLlC5cqVUy2fUr16ddV+Jycnnjx5wrNnz6hWrVqqMqpXr46vry9Pnz7F1taWK1euAKQ61sjIiMqVK6v2v69q1ars3r1bawtAitQkiYk0dezYkdGjRxMeHo65uTmbN2/G0dFRtR7Qf+3fvz9Fk2C/fv0YNmwYS5cuZdy4cSmaJZ8/f87s2bMZNGhQutf39fWlV69etGjRgiVLlqia/8aPH69amyj5en379qVPnz7MmjWLXr16YW1tTePGjZk9ezY1atTA29tbrd/5zp077Ny5U3VX0b17dypUqMDatWtVSWz16tUEBQWlaArs27evKtln5tatW0DmCTXZtGnTePv2LUeOHFGtKderVy/c3NyYMGECXl5eGBgY4OvrS7ly5Vi3bl2a5VSoUAEXFxe2b9+Ol5cXDg4Oqn379u3j1atXbN++PcWabBMnTsw0vtu3b5M/f34sLCzU+n0gaYmcvXv3qu7oLSwsmDhxIjdv3qRcuXJql5ORp0+fpjn4yM7ODoAnT56ojoN/F258n62treoYW1vbTI+9ePFiqu3Fixfn9evXPH36VHVtoV3SnCjS1Lp1axQKBbt27SI+Pp6dO3emaJr5r+SEkpCQQHh4OKGhodSpU4eoqKgUTTeQ9Jdr79690y1r165ddO/enQ4dOrBs2TJVAgsPD+fEiRO0adOG6OhoQkNDVf81adKEyMhILl26lOXf2dHRMUWzmKmpKdWrV0/R1HXkyBGsra1p06aNapuhoSH9+/dX6xqRkZEAmJmZZXpsQkICR48epVmzZikWRc2fPz9ffvklISEhXL9+XVXeo0eP0vwizUxyLAcPHiQuLk6jc8PCwjR+XKFPnz4pmqTd3d0B1G5SVEd0dHSaK3knb0tumo6Ojk6x/X3JfXrJx2R2bPL+9yUn99DQUI1/B6EeSWIiTQUKFMDDw4M//viD48eP8+LFiwyT2JkzZ2jevDl2dnYUL14cR0dHBg4cCJBqNJytrW26nf4hISF8+eWXeHh48Msvv6T4srt79y5KpZI5c+bg6OiY4r+hQ4cCfNCAk7RWNTY3N+fVq1eq1w8fPqREiRKp+gUdHR3VukZywkhOZhl5+fIlUVFRaa7qXaZMGQAePHgAwIgRI8iXLx+NGzemcuXKjBw5MtUovPTUrVuXli1bMmfOHEqWLIm3tzdr1qxJNWw8Pe+PylPHf+s5OQm+X88fKnfu3Lx79y7V9uRtye+/3Llzp9j+vuREl3xMZscm739fct1IU2L2keZEka6OHTvSq1cvIKkfIL0v6vv379O2bVtKlizJrFmzKFq0KKamply5coVvv/021eCOtD7syaysrLC3t+f48eOcPXtW1Z8FqMpJXpU4LR+yyN5/+0+SqfMlre4XeZkyZdi7dy83btzAzc1No/gyul7ZsmU5f/48R44c4ejRo/j6+rJ69Wr69OnDjz/+mGFZCoWC9evXc/HiRQ4ePMiJEycYMWIE8+fP5+jRo1hbW6d7bsGCBTUesv8h9awuGxsb1Wro70tuRkxu2ktuGkzr2OTmw+RmxfeP/e8Ak+Qmx/9KfiYvuQ9SaJ/ciYl0eXp6kj9/fvz9/TO8C9u/fz8xMTH8/vvv9O3bF09PTxo0aJClWTFMTU35/fffKV++PN7e3imaB5O/OIyMjGjQoEGa/yV/4WbXX77FihXj3r17qRJz8tLsmWnevDkAv//+e6bHFipUiLx583L79u1U+5KbaO3t7VXb8ubNS+vWrVm4cCFXr16lY8eOrF69msePHwOZ10m1atWYNGkShw8fZsuWLTx8+DDdPrZkZcqUISIigrCwsEx/n48peaBFQkJCiu0XLlwAUD1nVrhwYaytrdNshr1w4QI2Njaq5JQ8eOW/x8bHx3PlypU0n127d+8eBQoUSLMfTWiHJDGRLlNTU+bPn8+4cePo0KFDuscl/2X9/l/S7969Y9myZVm6br58+diyZQsODg60a9eOGzduAEl3afXq1WPNmjWEhISkOu/9psTkPjptz07RpEkTnj9/zs6dO1XbEhISWL58uVrnV69eHQ8PDzZs2MCuXbtS7U9ISGDhwoU8evQIQ0NDGjdujK+vL3fu3FEdExkZyerVqylatKhqgMh/k4iRkZFqX3IdpFcn4eHhqe6Ckr+QM6u/WrVqAWTbLCNZ1bp1a8LCwti6datq27t371i7di3Ozs44OTmlOPbIkSMpZui4e/cux44dSzFgp3Tp0jg7O7N27VpiY2NV2zdv3kx4eHiag3sCAgKoWbOmNCdmI2lOFBnKKHkla9y4MSYmJnTu3JnevXsTGxvL77//nubzZOoyNzdnx44deHl50bZtW/bv34+joyMLFizA09MTd3d3evXqhaOjI6GhoVy5coVjx47x8OFDIKmPKn/+/KxatYp8+fKRL18+ypUr90HNjQC9e/dm+fLlDB48mICAANUQ++RZGdT5svLx8aFDhw706tULDw8PGjRogJmZGcHBwezatYu7d++q6v2bb77hxIkTNG/enH79+qmG2IeEhLBmzRpVHbdt2xYrKytcXV2xtrbm3r17LFu2DGdnZ8qWLQugGnk4ffp02rdvj4mJCfXq1WPLli2sWLGCFi1aUKJECaKjo9m4cSOGhoaZjrqsWbMmVlZWHD9+XDWCMzv9/vvvqn9jSBoq/8MPPwBJz2sl35m2atUKd3d3hg8fTmBgIEWLFmXTpk3cu3ePbdu2pSjz66+/ZteuXbRs2VL1qISPjw8WFhaMGjUqxbHff/897du3p1WrVnTu3JmQkBAWLVpEnTp1VI87JHv69Cn//PMPAwYM0Ho9iH9JEhMfrFSpUmzcuJFp06bx7bffYmlpSefOnalTpw5t27bNcrmFChVi586dfPHFF7Ru3Zr9+/dTqlQpTpw4wdy5c9myZQsvX77E0tKSMmXKMH36dNW5pqamLF26lOnTpzN69Gji4uIYN27cByexvHnzsmfPHsaNG6d6dq1ly5ZMmjQJT09PtWapsLS05ODBg6xevZpt27Yxe/ZsoqOjsbOzUz2gXLhwYQCcnJw4ePAgU6dOZfHixcTGxuLi4sLvv/+eol+wT58+bNmyBR8fHyIjI7G1taVbt26MGTNGlehq1KjB5MmTWbNmDUOHDiUxMZE9e/bg7u7OpUuX2LFjB8+fP8fMzIyKFSsyd+7cTOfINDY2xtvbmx07djBt2rRsv+NYv349/v7+qtcBAQEEBAQA4OrqqkpiCoWCTZs2MX36dNauXUtkZCTOzs5s3rw51YPZtra27N+/n0mTJqlm53Bzc2PmzJmp+rkaNmzIH3/8wcyZMxk/fjxmZmb06NGDKVOmpPrdd+/ejamp6Qd9BkTmFOHh4drrTRXiM7Vnzx569OjBwYMHcXV11XU4H9XDhw+pXr06a9eupVmzZroORy8olUrc3d2pV69eimmrhPZJn5gQGvrv80AJCQksWbKE/PnzazRzxaeiWLFi9OnTR2Zsf4+vry/BwcEZTqkmtEPuxITQUPv27bG2tqZKlSrExMSwa9cuLl68yNSpUxk+fLiuwxPisyJJTAgN+fj4sH79eh48eEBcXByOjo7079+fPn366Do0IT47ksSEEELkWNInJoQQIseSJCaEECLHkiQmhBAix5IkJoQQIseSJCaEECLHkiQmhBAix/o//qGkOZBwsYEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "toy_panel = pd.DataFrame({\n",
    "    \"mkt_costs\":[5,4,3.5,3, 10,9.5,9,8, 4,3,2,1, 8,7,6,4],\n",
    "    \"purchase\":[12,9,7.5,7, 9,7,6.5,5, 15,14.5,14,13, 11,9.5,8,5],\n",
    "    \"city\":[\"C0\",\"C0\",\"C0\",\"C0\", \"C2\",\"C2\",\"C2\",\"C2\", \"C1\",\"C1\",\"C1\",\"C1\", \"C3\",\"C3\",\"C3\",\"C3\"]\n",
    "})\n",
    "\n",
    "m = smf.ols(\"purchase ~ mkt_costs\", data=toy_panel).fit()\n",
    "\n",
    "plt.scatter(toy_panel.mkt_costs, toy_panel.purchase)\n",
    "plt.plot(toy_panel.mkt_costs, m.fittedvalues, c=\"C5\", label=\"Regression Line\")\n",
    "plt.xlabel(\"Marketing Costs (in 1000)\")\n",
    "plt.ylabel(\"In-app Purchase (in 1000)\")\n",
    "plt.title(\"Simple OLS Model\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ed02cca",
   "metadata": {},
   "source": [
    "对因果推理了解很多后，您决定运行一个固定效应模型，将城市标识作为虚拟变量添加到您的模型中。 固定效应模型控制了在时间上保持不变的城市特定特征，因此如果一个城市对您的产品不太开放，它将捕捉到这一点。 当您运行该模型时，您最终可以看到更多的营销成本会导致更高的应用内购买。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b34af3b8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAE0CAYAAACirQ3aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABpEUlEQVR4nO3dd1QU59fA8e9SRQFBRFABC2LB3hEVG1JEY2yxx1hiNCZRExU1xaiJWGM0xRoTW9RYYhfsit3Yu9jFLkrv7L5/8GNfcQAXBRb0fs7JOWHu7MzdAffuzNx5HlV4eLgGIYQQogAy0HcCQgghxOuSIiaEEKLAkiImhBCiwJIiJoQQosCSIiaEEKLAkiImhBCiwJIiJvTCysoKPz8/faehFRAQgJWVFcuXL8/V/Zw6dYoOHTpQoUIFrKysqF69ujZ2/fp1evXqReXKlbG2tsbKyipXc3kbDR48GCsrK4KDg99oO2l/D2+6HZH7jPSdgHi7vOqDd9KkSXz66ad5k0wuWr58OUOGDMlyHUtLS+7cuaP9OTIykg8++IDIyEi6du1KyZIlKVq0KAApKSn07NmTK1eu0KVLF8qWLYtKpcrV95DGysoKR0dHzp07l63XBQQEMGXKFAC6dOnCggULMlzvyJEj+Pj4AFCiRAmuXr36ZgkL8QIpYiJX+Pv7Z7i8fv36ABw7dgwzM7O8TClXVKtWLdMzykKFCqX7+eTJkzx58oS+ffsyc+bMdLHbt29z+fJlWrVqxfz583Mt39xgZGTExo0bmTp1KtbW1or44sWLMTIyIjk5WQ/ZibedFDGRK8aMGZNlvGLFinmUSe6qXr36K99rmgcPHgCpZyPZieV33t7ebNmyhZUrVzJ48OB0sYiICDZs2ICPjw+bN2/WU4bibSb3xIRevHxP7M6dO5QtWxYnJydu3bqVbt3ExESaN2+OlZUVmzZtShfbv38/3bp1w9nZGVtbW6pVq8ZXX33Fo0ePMtzv6dOn6dSpEw4ODjg6OtK+fXuOHj2a4+/vRbdv38bKykr7AT9lyhSsrKy09+BePBYrVqzQxgICArTbUKvVLFmyBG9vb5ycnLCzs6NRo0b89NNPJCYmZrjf69evM3ToUGrWrImdnR3lypWjZcuWTJs2DYDg4GDt5d+7d+9q9/tirrpo0aIFjo6OLFmyRBFbvXo1sbGx9OnTJ9PXazQa/vrrL1q1aoWDgwMlS5akSZMm/PLLL5m+t7179+Lr60upUqUoW7YsPXr04MqVK1nmeePGDT7//HOqVatGiRIlcHZ2pmfPnpw+fVrn9yryHzkTE/mCk5MTv/32Gz179qRv374EBQVhYmICwDfffMPp06cZOHAg7dq1077m559/5vvvv8fa2hovLy/s7Oy4cOECf/zxB9u2bWPHjh2ULl1au/7Ro0d5//33SUhIoF27djg7O3PhwgXatWuHh4dHrr23okWL4u/vz7lz59i6dSuNGzemSZMmQOqZnL+/P3fu3GHFihXpLk+mrZOcnEyvXr0IDAykQoUKdOrUCVNTUw4ePMiECRPYt28fa9euxcjo//8579y5kw8//JC4uDiaN29Ohw4diImJ4dKlSwQEBDBy5EicnJzw9/dnypQpWFpapitcLzacvIqBgQG9e/dm0qRJHD9+XHvJGFIvJTo4ONCqVatMXz9w4EBWr15NqVKl6NGjB8bGxgQGBvLtt9+yc+dOxXvbsGEDffv2xdjYmPfff59SpUpx5MgRWrduTbVq1TLcx759++jZsyfx8fF4e3vj7OzMgwcP2LRpEzt37uTvv//OMkeRf0kRE7nixbOINHZ2dvTr1y/T1/j5+TFo0CDmzp3Ld999x+TJk9m0aRPz58+nZs2a/PDDD9p1Dx48yPjx46lfvz6rV69O11CycuVKBg0ahL+/P8uWLQNSv+1/9tlnxMXFsXjxYtq3b69df8GCBYwcOfK13ue5c+cyfK+QWgjatm2LlZUVY8aMYfny5WzdupUmTZqkuwRZo0YNgoODWbFiRYaXJ2fOnElgYCAff/wxkydPxtDQEEg9Oxs+fDiLFy9m4cKFDBo0CICwsDD69etHfHw8q1evxtPTM932QkNDAShTpgxjxoxhypQpFC1aVOfLohnp1asXU6ZMYfHixdoidurUKc6dO8fo0aMxMMj4os+aNWtYvXo1VatWZdu2bVhaWgIwbtw4OnfuzL59+/j999/54osvAIiOjmbYsGGoVCq2bNlCvXr1tNv69ttv+eWXXxT7iIiI0Ba9nTt3UrlyZW3sypUrtGrViiFDhnDmzBlMTU1f+xgI/ZAiJnJFWtfai6pVq5ZlEQOYMGECR48eZe7cuZQtW5aAgAAsLS35888/tWdmAHPnzkWj0TBz5kxFR2S3bt34/fff2bZtG5GRkVhaWnL06FFCQkJo2LBhugIG0L9/f+bMmcONGzey/T7Pnz/P+fPnM4x1796dtm3bZnubL1Kr1cydOxdbW1sCAgK0BQxSz4AmTJjAkiVLWLVqlbaI/f3330RGRtK3b19FAQNwcHB4o5wyUqpUKTw9Pfn3338JCAjAwsKCxYsXY2BgQK9evTJ9XdqXjHHjxmkLGICJiQmTJk2iSZMmLF68WFvEtm7dyvPnz+nSpUu6AgYwatQoFi9eTGRkZLrlK1eu5NmzZ0yePDldAQOoVKkSH374Ib///jt79+7F29v7jY6DyHtSxESuCA8Pf63XmZiY8Oeff+Lh4cHo0aMBWLRoEeXLl0+33tGjRzEyMmLTpk2K+2SQeh8tJSWFGzduUKtWLc6cOQNA48aNFesaGBjg5ub2WkWse/fuzJkzJ9uv09W1a9cICwujXLly2ntZLzMzMyMkJET783///QeAl5dXruWVkT59+hAUFMSaNWvo0qULa9euxdPTM8uimfZ7adq0qSJWrVo1bG1tuX79OtHR0Zibm2f5e7SwsKBGjRocOHAg3fK0e54XLlzI8Kz52rVrAFy9elWKWAEkRUzkO2XLlqVmzZoEBwfj4OCQ7j5YmmfPnpGcnJzhGd+LoqOjAbTfzm1tbTNcL792BT579gyAmzdvvvK9pomIiABSz47ykre3N6VKlWLJkiUYGhoSFRXFhx9+mOVr0s6UM3vcws7OjidPnhAZGYm5uflr/R7TjuHSpUuzzCUmJibLuMifpIiJfOfXX38lODgYGxsbQkNDGT9+fLr7YZD6IHFSUhJ3797VaZtpl6qePHmSYfzx48dvlnQuScvbx8eHlStX6vSatAeoHzx4QM2aNXMtt5cZGhrSo0cPpk+fzuPHj7G3t9c+5JwZS0tLnj9/TlxcXIaFLK3LNO04vM7vMe01e/fupVatWjq/H1EwSIu9yFdOnDjB+PHjcXJy4vDhw1SvXp3ffvuNoKCgdOvVr1+fqKgonUeZSPswP3jwoCKmVqs5cuTImyefCypWrEjRokU5ceJEpu3mL0trrNi+fbtO6xsYGKBWq187xxf17t0blUrFvXv36NmzZ7quwoyk/V5evgQIcPHiRZ48eUKFChUwNzdPt35Gv8eoqCjOnj2rWJ52PA4fPpy9NyMKBCliIt8IDw+nb9++QOp9sBIlSvDXX39hbm7O4MGDuX//vnbdtCGfhg0bxr179xTbio+PT/eh1bBhQ1xcXDh69CgbNmxIt+4ff/zxWvfD8oKRkRGDBg3iyZMnjBgxgtjYWMU6YWFh6T68e/TogaWlJYsXL2bPnj2K9V8+XjY2Njx9+pS4uLg3zrdMmTL8+++/LFu2TKdnzXr37g2kNvSkXfoFSEpK4uuvvwZId0myTZs2WFlZsW7dOu29vzRTp05VNHVAaueklZUV06ZN49ixY4q4RqPh8OHDOn9JEPmLXE4U+cbnn3/OnTt3mDhxorbzzNnZmZkzZzJgwAD69+/P5s2bMTQ0xMPDg4kTJzJu3Djq1q1L69atKVu2LPHx8dy9e5dDhw7h5OSk/YavUqn45Zdf6NChA3379k33nNiePXvw9PRk586d2c45qxZ7gOHDhyuGn8qukSNHcvHiRZYsWcL27dvx8PCgdOnSPH36lJs3b3LkyBEGDBhAjRo1AChWrBiLFi3iww8/pGPHjrRo0YKaNWsSExPD1atXCQ4OJiwsTLv9Fi1a8M8//9CpUyfc3d0xNTWlWrVq+Pr6vla+zZs313ndTp06ERgYyOrVq3Fzc8PPz0/7nNi1a9do1qxZumJobm7OrFmz6Nu3L35+fnTo0IFSpUpx+PBhLl68iLu7O4cOHUq3D2tra5YsWUKvXr3w8vLCw8ODypUrY2xszL179/jvv/8IDQ3l1q1b6TpgRcEgRUzkC/PmzWPTpk14eXnx2WefpYt17tyZ/fv3s2TJEgICAvjmm2+A1KLn5ubG3LlzOXz4MIGBgZibm1OyZEm6dOlCx44d023Hzc2Nbdu2MXHiRHbt2sWuXbuoW7cumzdvZteuXa9VxLJqsYfUUdXftIgZGRmxZMkS1q5dy/Lly9mxYwfR0dEUK1YMR0dHhg8fTrdu3dK9xtPTk7179/Lzzz+zb98+goODsbCwoHz58owdOzbdupMnT8bAwIA9e/Zw9OhRUlJS6N69+2sXseyaN28e7u7uLF26lKVLl6JWq3F2dmbChAkMGjQIY2PjdOu3b9+etWvXMmXKFDZs2ICJiQnu7u7s2LGDmTNnKooYgIeHBwcPHuTXX39l165dHDt2DCMjI+zs7Khfv76ixV8UHKrw8HCNvpMQQgghXofcExNCCFFgSRETQghRYEkRE0IIUWBJERNCCFFgSRETQghRYEkRE0IIUWBJERNCCFFgSRErAF6cZkO8mhyv7JHjlT1yvLInt4+XFDEhhBAFlhQxIYQQBZYUMSGEEAWWFDEhhBAFlhQxIYQQBZYUMfFW0WhkUgYh3iUyn5h4K6hT1Jycf5yIOxE49Sun73SEEHlE5yIWFhbGkSNHuHr1KmFhYahUKmxsbKhYsSINGzbExsYmN/MUIlNRD6LYPmwroUfuAmBc3gQXFxc9ZyWEyAtZFrGEhARWr17N8uXLOXr0aKaXalQqFQ0aNKBnz5588MEHmJqa5kqyQrwsZOtVdo3eTkJEvHbZuZmnqdWmDpalZaZeId52md4T+/PPP6lduzZffvkllpaW/PDDD2zbto1Lly7x8OFDHjx4wKVLl9i6dSsTJkzAwsKCr776itq1a/Pnn3/qtPODBw/SrVs3qlSpgpWVFcuXL1esc+3aNXr16oWTkxMlS5bEw8ODK1euvP47Fm+FxJhEdo4KYuvgjekKGIDK0IDIuxF6ykwIkZcyPRObNm0agwcPpnfv3lhZWWW4jr29Pfb29jRq1IjPPvuM8PBwli5dyvTp0+nbt+8rdx4TE4Orqyvdu3dn0KBBivitW7fw9vamW7dubNy4ESsrK65evUqRIkV0f4firfPwzAOChm4l/OZzRcyxSRlcPquEg5ujHjITQuQ1VXh4eIbXCJOSkjA2Nn6tjb7Oa0uXLs3UqVPp2bOndtmAAQNQqVQsWLDgtfJ4W4SEhMg9HlKbN07MPc6Rnw6iTlanixkYG9B4VFNqD6jHtevX5Hhlg/x9ZY8cr+zJ7eOV6eXE1y1gb/raNGq1msDAQCpVqkSnTp1wdnamRYsWrFu37o23LQqeqPuRrOvxD4emBisKmLVzMbqu70mdgfVRGaj0lKEQQh907k6MiYnh0qVLPHjwgPj4eAoVKkTJkiWpUqVKrlzee/LkCdHR0fz000+MHTuWcePGsX//fj7++GMKFy6Mj49Phq97W0eYflvfly4e7LvHuZlnSI5OUsSc2palyqCqRJhGEhESqV3+Lh+v1yHHK3vkeGXPmxyvV53FvbKIXb16lQkTJrBjxw6SkpLSdSiqVCqMjY3x9PTk22+/pXLlyq+d6MvU6tRv223atOGzzz4DoEaNGpw+fZqFCxdmWsTextP8d/XyRWJMIvvG7ebi6vOKWCFrMzyneuPsVUERe1eP1+uS45U9cryyJ7ePV5ZF7PTp07Rr1w4TExP69OlDvXr1sLe3p1ChQsTHx/Pw4UOOHz/OunXraN26NZs2baJWrVo5kpiNjQ1GRkZUqlQp3fKKFSvKJcV3wMPTDwgcuoWIW+GKmFPTMnjN8KWInXneJyaEyFeyLGLfffcdjo6ObNmyBWtr6wzX+eCDDxgzZgx+fn6MGzeODRs25EhiJiYm1KlTR3Eaeu3aNRwdpfPsbaVOUfPfnGMcnXlIce/L0MQQd/+m1O5XV+59CSGAVxSxkydPMn78+EwLWJpixYrRr18/xo8fn62dR0dHc+PGDSD18mFoaChnz57F2toaR0dHvvjiC/r27Yu7uzseHh4EBwezbt26DJ8nEwVf5L1Itg/fyr2joYpYsQrF8JndFtuqJfSQmRAiv8pyAGATExMiIyOzWkUrKioKExOTbO381KlTeHh44OHhQVxcHAEBAXh4eDBp0iQA2rZty88//8wvv/yCu7s78+bNY+7cuXh7e2drPyL/u7r5Mn/7LM6wgNXoXYtum3tLARNCKGR5Jubp6cmvv/6Ku7s7DRs2zHS9o0eP8uuvv9K6dets7bxp06aEh4dnuU7Pnj3TPTsm3i6J0YnsHbeLS2suKGJmxczwnOZDeU9nPWQmhCgIsixiP/74I2fPnsXX15dq1apRt25d7OzsMDU1JSEhgUePHnHixAnOnz+Pi4sLP/zwQ17lLd4CD0/9r3njdrgi5uRRFq/pPtK8IYTIUpZFzNbWlr179/LHH3+wfv16Vq5cSXz8/49TZ2pqSvXq1ZkwYQL9+/fHzMws1xMWBZ86Rc1/vx/jyMyDaFLSDxhjaGJI49Ee1OpbR5o3hBCv9MrnxAoVKsSQIUMYMmQIGo2G58+fExcXh5mZGdbW1qhU8kEjdBcZGkHQ8G3cP5ZB84aLDT6/tMW2iq0eMhNCFETZmhRTpVJRrFix3MpFvOWubLzM7q93kBiZoIjV7FObJmM9MCr05kOWCSHeHToVsaSkJA4cOMCZM2cUw07VrFmTJk2a5Mh4ieLtlBCVwL5xu7i09qIiZmZjRutpPpRrJc0bQojse2URW7VqFd999x1PnjzRDjllbGxMUlLqOHYqlYrixYszfvx4unfvnrvZigLnwcn7BA3dQsQd5fxeZZqXo/V0H4rYytQ6QojXk2URW7NmDYMGDaJJkyZMmzaNevXqUbJkSVQqFRqNhgcPHnD8+HEWLFjAkCFDMDY2pnPnznmVu8jH1Clqjv96lKOzDimbN0wNaTKmGTU/qi33VIUQbyTT+cQAGjdujIODA6tWrXrlhrp27UpoaCgHDx7M0QRFwRtwNPJuBEHDt3L/+D1FzKZiavNG8cq517xR0I6Xvsnxyh45Xtmjt/nEAK5fv46fn59OG/Lz8+P69es5kpQouK5suMRy38UZFrCaH9Wm26ZeuVfA1CkY/rePkrvXYnDpFGgy/X4mhHhLZHk5sXjx4ly6dEmnDV28eJHixYvnSFKi4EmISmDvt7u4/G8GzRvFC6c2b7Qsn3sJxMVgNm0kBndCMEtKQnNyHymVahI/9AcwzFYTrhCiAMnyTKx79+4sWLCAGTNmZDo8VHh4ODNmzGDhwoX06NEjN3IU+dz9/+7xt++SDAtY2Rbl6BnYJ3cLGGCyYg6G1y+iSms4io/D8NxxjHblzKwKQoj8KcuvqP7+/ty7d48ffviBgIAAypUrpxh26ubNm6SkpNC1a1dGjRqVV3mLfECdrObYr0c4Nvtwhs0bTcc2o0afvGneMAy9oVimUqdgdP44yV6dcn3/Qgj9yLKIGRkZ8fvvvzNw4EDWr1/P2bNnefjwoXbEDgcHB9q2bUv79u1zbDJMUTBE3o0gcOgWHpy4r4jZVC6Oz2w/ilfKu5E3NMYZz6CgMS2UZzkIIfKeTjcLatWqJUVKaF3+9yJ7vt1JYlSiIlarXx0a+3tgVChv70Mlu7fG8PZVVHGx2mVqS2uSfLvlaR5CiLwld7yFzhIiE9jz7U6urFc2+xS2LUzr6b6UbV5OD5lBcjM/VJHhGB3ZSUpEOIY2tiR5dUZdvrJe8hFC5A2diti5c+dYv359hsNO1apVi/fee48aNWrkdq5Cj+4fv0fQsC1EhionSS3bsjytp3lTuLh+R95IateTJL/uXLt8kQpVqoI8SC3EWy/LIpacnMzw4cNZvnw5AGXLlsXe3h47Ozvi4+O5ceMGu3fv5qeffqJ79+7MmjULIyM5uXubqJPVHPvlMMdmH0Gjfrl5w4imXzejxoe18s/IGwYGaIxNpYAJ8Y7IsuJMmzaNv//+m1GjRvHxxx9jY2OjWOfZs2fMmzeP6dOn4+DgwJgxY3ItWZG3Iu6EEzR0Kw9OKps3ile2xeen1ti4lpSCIYTQmyyL2N9//83HH3/M6NGjM12nWLFijBkzhufPn7N8+XIpYm8BjUbD5X8vsvfbXSRGK5s36npa0KJEMMYLN6OxsCKpVQeSm7XRQ6ZCiHddlg87P3nyBFdXV502VK1aNZ4+fZojSQn9SYiIJ/CLLWwfvk1RwArbFqHj2Ap4GW/E9P41DJ49wfB2CCar52Fw5YyeMhZCvMuyLGLOzs5s27ZNpw1t3boVZ2eZE6ogu3c8lOVtlnB142VFrJynMz2D+uAScQCDuJh0MYOoCEwCV+dVmkIIoZVlERs6dCiBgYF07NiRbdu28ejRo3TxR48esXXrVjp27Mj27dsZOnRoriYrcoc6Wc3hGQdY+8Eqol7qPjQ0NaLFD560W/g+hW0Ko0pSzsoMQGImy4UQIhdleU/sgw8+IDk5mXHjxtGzZ0/t8hcnxdRoNNjY2DB79mw++OCD3M1W5Ljw2+EEDd3Cw1MPFLHirrb4zPLDpuL/D+ycUsYFw6vn0q2nQUVK5Zq5nqsQQrzslf3wPXr0oHPnzgQHByueE7O3t6dWrVo0adIEU1PTvMhX5BCNRsPldakjbyTFJCnitQfUxX1UU4xM0/+JJHb5GMNbIRjcuoIqKRGNaSFSKlQlybdrXqUuhBBaOj3UZWJiQqtWrWjVqlVu5yPyQEJEPLu/3sHVTVcUscK2RfD6yZcyHmUzfrGpGXFjZ2F46hCG186T7FoXdbV60mYvhNALeTL5HXPv6F2Chm8l6l6UIla+tTOeU70xK1Y4640YGJBStwkpdZvkUpZCCKGbLBs7siMwMJAhQ4Zk6zUHDx6kW7duVKlSBSsrK+3IIBkZOnQoVlZW/PLLL2+a6jspJSmFQ9MOsLbbP4oCZlTIiJaTWtN2wfuvLmBCCJGP5FgRO3/+PCtWrMjWa2JiYnB1dWXy5MmYmZllut6GDRs4efIkJUuWfNM0C5boCEz+/o1yq37B5O/fIFo5bqEuwm89Z3XnFRz/VTl0lK1rCbpv7k31njXzz9BRQgihI71eTvTy8sLLywuATz/9NMN17ty5w+jRo1m/fj2dO3fOy/T0ShXxjEKTh2N4/zYmACFnMTx3jPjRM9EULabTNjQaDZfWXGDvuF0ZNm/U+bgejUY2UTRvCCFEQZHlp1d2znxSUlLeOJmXJScnM2DAAEaMGEGlSpVyfPv5mcnqBRjev51umeH92xivWUhi/1fPoB0fEc/usTsI2axs3ihSogheM9vg1KRMjuUrhBD6kGURS0pKokyZMjRs2PCVG7pw4QLnzp175XrZERAQgLW1Nf3798/R7RYEqifK57YADDJZ/qLQI6nNG9H3M2je8KqA5xQvufclhHgrZFnEKleuTOHChfn9999fuaHp06fnaBE7cOAAf//9N8HBwdl6XUhISI7loE9lMcA6g+VRGHArk/eoTlYTsuQy11eEQPpbXxiYGuI6uBqOfmUIDbsHYTmfc37ytvwd5BU5Xtkjxyt73uR4ubi4ZBnPsojVrl2bNWvWkJycnOfzhAUHB/Pw4cN0lxFTUlIYN24cc+bM4eLFixm+7lVvuKBQ9RyCeuZoDJ490S5TFyuBSa8huDiUV6z//OZzgr7awqMzDxWxEtXs8J7VhmIVlFPpvI1CQkLemr+DvCDHK3vkeGVPbh+vLCtTx44dSUlJISwsDDs7uyw35OvrS6lSpXIssQEDBtC+fft0yzp16kSnTp3o06dPju0nv9I4ORP/2XhM1i8m4ckjTG3tSOzwEZqXCphGo+HiP+fZ9/1ukmKVzRt1P6lPoxFNMDQxzKvUhRAiz2RZxFq0aEGLFi102lDVqlWpWrVqtnYeHR3NjRs3AFCr1YSGhnL27Fmsra1xdHTE1tY2fbJGRtjZ2b0z34LUzq7EfzUl028y8eFx7Bqzg2tbrypiRezM8frJV5o3hBBvtRx7Tux1nDp1Cg8PDzw8PIiLiyMgIAAPDw8mTZqkz7QKhNDDd1jusyTDAubs7ULPoD5SwIQQbz29PiDUtGlTwsPDdV4/p7sfC6KUxBSOzDzIf3OOKZo3jMyMaDauJVW7VZcHl4UQ7wR5yrUAeX7jGYFDt/D47CNFrER1O3xm+WHtrNuD0EII8TaQIlYAaDQazq88x77vd5Ecl5w+qIK6gxrQ6MvG0rwhhHjnSBHL5+LD4zg14TgPg5UPOZvbm+M1sw2O7k56yEwIIfRPilg+dvfQHbYP30r0w2hFzNnHhVaTvTCzznzgZCGEeNtJEcuHUhJTODzjACfmHc+4eeP7VlTtWk2aN4QQ77xsFbGQkBCWLVvGrVu3eP78ORpN+k9YlUrFxo0bczTBd83z688I/GILj89n0LxR43/NG+WleUMIISAbRWzt2rV88sknGBoa4uLigpWVlWKdl4ua0J1Go+HCynPsG787w+aNeoMb4DZcmjeEEOJFOhexSZMm4erqytq1axUjaYg3E/c8jl3+27kepBwk07ykBVW/qoFbl0Z6yEwIIfI3nYvYvXv3+OGHH6SA5bA7B26z/cttxDxSNm9UaFORVgGtufskVA+ZCSFE/qdzEatYsSJhYW/n/B0GNy5hsnYRqqgINOYWJL7XG3XlWrm6z5TEFA5PP8CJ+crmDePCxjQb3xLXLv9r3niS8TaEEOJdp3MR++677/jss8/o1KkTFSpUyM2c8pTq/m0K/ToOg7DH2mUG924T//kE1BVcc2Wfz66FEfjFFp5ceKyI2dW0x3uWH9blMppNTAghxIt0LmLbtm3D1tYWd3d3PDw8cHBwwNAwfZOBSqVi+vTpOZ5kbjLZsDhdAQMwCH+KyaZlxA/P2YGINRoN5/8+y/4Je0iOVzZv1P+0IQ2Hu2NoLM0bQgihC52L2KJFi7T/v2vXrgzXKYhFTBUdmfHyWOU9qjcR9yyWnf7bubH9miJmXsoC75ltcHBzzNF9CiHE207nIvb8+fPczENv1PaOcP4/5fLiWU8Cmh23g2+x48ttxDyOUcRc2lai5aTWFCpaKMf2J4QQ74p3fsSOxA59MbxyFsO717XLUkqVIfGDT95428kJyRyedoCTC5RF0riIMc0ntKJKp6oy8oYQQrymd76IYW5J3NezMd6yAsN7t1DblSaxbU8wt3yjzT4L+V/zxsUMmjdq2eMzyw+rstK8IYQQbyLTIlajRg0MDAw4fvw4xsbG1KhR45VnDCqVitOnT+d0jrnPrAhJnQeQlAOb0mg0nFt+huCJexXNGyoDFfWHNKTB0EbSvCGEEDkg0yLWuHFjVCoVBgYG6X4WmYsNi2WXfxA3dlxXxCxKW+D9sx+lGzjoITMhhHg7ZVrE5syZk+XPIr3b+2+x/cttxD5RNm9UbFeJlj+2xlSaN4QQIkfJPbE3lJyQzKEpwZz644QiZlzEmBYTPanc0VXOYoUQIhdkWsSuXbv22iNzvMlrC5Kwq08J/GILTy8px4Wyr10S71l+WJWxyvvEhBDiHWGQWcDNzY3+/ftz8OBBnTak0WjYv38/ffr0oVGjt3vEdY1Gw5klp1jRdpmigKkMVDT4wo3Oq7tJARNCiFyW6ZnY7t27mThxIm3btsXe3p6mTZtSu3ZtypQpg5WVFRqNhvDwcG7fvs3p06fZv38/jx8/plWrVpmO6PE2iA2LZefIQG7uuqGIWThY4v1zG0rXl+YNIYTIC1m22K9evZpLly6xbNkytmzZwurVqwG093fSJsEsU6YMHTt2pFevXri65s6gufnB7X032f7VNmKfxCpiFd+rTMsfPKV5Qwgh8tArGzuqVKnCjz/+yI8//sjDhw+5evUqz549A6BYsWJUqlQJO7ucG6IpP0qOT+bg1GBOZ9C8YWJuom3eEEIIkbey1Z1ob2+Pvb19buWSL2XVvFGyTim8Z7WhqJNV3icmhBBCWuwzo9FoOLv4FMGT9pOSoBx5o8EXbjT4vBEGRpn2xgghhMhlev0EPnjwIN26daNKlSpYWVmxfPlybSwpKYlx48bh7u5OqVKlqFSpEgMGDODu3bu5nlfs0xg29l3H3nG7FQXM0sGSzv90w214YylgQmfJag0RiWrtfWQhRM7Q66dwTEwMrq6uTJ48GTMzs3Sx2NhYzpw5w4gRI9i3bx9///039+7do3PnziQnJ2eyxTd3a88Nlnkv5taem4pYpfer0GNbH0rVL51r+xdvF41Gw8QTEbivf0zDdY9otvEJq68rG4OEEK9Hr5cTvby88PLyAuDTTz9NFytatCjr169Pt2zmzJm4ublx5coVqlatmqO5JMcncyBgH2f+OqWImVj8r3mjgzRviOz5/WI0cy5EE5uS+vPDODVfH4vA1dqIqsVM9JucEG+BAnU9LCoqCgArK6sc3/bl9ZcyLGAl65aix9YPpYCJ17LxVry2gKV5HK9m9rmcnTlciHdVgWnsSExM5JtvvsHHx4fSpTO/nBcSEvJa2zepbUrxurY8PZHahagyUFGhdyWce7jwOOEJj0OU3Yl56XXf17sqvxyvyFhTQDntzsPwKEJCwvI+oUzkl+NVUMjxyp43OV4uLi5ZxrNVxCIjI5k/fz779+/n6dOnzJ49m3r16vHs2TOWLl1K27ZtcXZ2fu1kM5OcnMzAgQOJiIhgxYoVWa77qjeclVJzSrHcezEm5ib4zPKjZN1Sr72tnBQSEvJG7+tdk5+OV60Hz7gUHZdumYkBdHMtjotLET1llV5+Ol4FgRyv7Mnt46VzEbt//z5t2rTh3r17ODs7c/XqVWJiUqcdKVasGEuWLOH+/ftMmTIlRxNMTk6mf//+XLx4kc2bN1OsWLEc3f6LzO3Mab+4I9bli2FqYZpr+xHvjoCGVoREJnM+LIl4NRQ1VtG8tCldKxTWd2pCvBV0LmLff/89kZGR7Nu3Dzs7O8Uo9X5+fmzfvj1Hk0tKSqJfv35cunSJzZs358nIIPY1S+b6PsS7w8rUgO1+tgTejedsWBJeDoWoYysNHULkFJ2L2M6dO/nkk09wdXXVDjv1orJly3L//v1s7Tw6OpobN1IH0lWr1YSGhnL27Fmsra0pWbIkffr04dSpU6xYsQKVSsWjR48AsLS0VLTkC5FfGahUtHEyo42T/M0KkdN07k6MjY3N8kwoNjYWtVqdrZ2fOnUKDw8PPDw8iIuLIyAgAA8PDyZNmsS9e/fYunUrDx48oHnz5lSqVEn737p167K1HyGEEG8nnc/EnJ2dOXHiBB999FGG8Z07d2Z7BPumTZsSHh6eaTyrmBBCCKHzmVifPn1YuXIlK1eu1J5xqVQqYmJi+Pbbb9m/fz/9+/fPtUSFEEKIl+l8JjZw4EAuXbrE4MGDsbCwAKBfv36Eh4eTkpLCJ598QteuXXMtUSGEEOJl2XpObObMmXTr1o1///2XGzduoFarKVeuHB07dqRRo0a5laMQQgiRoWyP2NGwYUMaNmyYG7kIIYQQ2aJzEYuLiyMmJobixYtrlz19+pQlS5YQHh5O+/btqVu3bq4kKYQQQmRE5yI2fPhwLl26xL59+4DUaVRatmypnd9rzpw5bNq0CTc3t9zJVAghhHiJzt2JR44cwdfXV/vzmjVruHv3LmvWrOHKlStUqlSJ6dOn50qSQgghREZ0LmKPHj1KN3r8tm3baNCgAa1ataJEiRL07NmTs2fP5kqSQgghREZ0LmJFihTRPnycnJzMoUOHaN68uTZuZmamne9LCCHeRo+eh3Lv2TXiE2V27vxC53titWvXZunSpXh4eLBt2zaio6Px8fHRxm/evEmJEiVyJUkhhNCn+MRY/t49iwdhd4hPiuXk3V3Ur9SCptX99J3aO0/nIvbNN9/QoUMHWrRogUaj4b333qN27dra+ObNm6X1XgjxVtpw6C9uPrys/fl59BMOXNhGRcda2FllPkmvyH06F7GaNWty/Phxjh49ioWFBU2bNtXGwsPDGTBgAI0bN86VJIUQQp8ePr+jWBYbH8XRizt4z/2jvE9IaGXrYWcbGxvatGmjWG5lZcXgwYNzLCkhhMhPVKgwuWyBxkhDUoVo7XIDA53bCkQuyfaIHQBRUVFERkZmOPWKo6PjGyclhBD5RVJcEoU3lUSzywZ1kSSef3oNjXky5mZWNK7q++oNiFyVrSK2ZMkSZs+erZ3IMiMZTZgphBAF0ZMLjwn8YjMx1xIBMIgxxnK9I8afxtO4qg/WFrZ6zlDofC68dOlShg4diqOjI9988w0ajYbBgwczfPhwSpQoQfXq1fnll19yM1chhMgTGrWGk/OPs7L9Mp5dS//F3OSGBV3KfU69Ss31k5xIR+ciNmfOHJo2bcq///6rnRjTy8uLb7/9liNHjhAeHk5kZGRu5SmEEHki+lE0//ZeQ/CP+1Anpb9lUrSsFY1mNaVUbelIzC90LmI3btygbdu2qS/6383MpKQkILWx48MPP2ThwoW5kKIQQuSN60EhLPdezN0DtxUx1w+q0WPrh1hVttZDZiIzOt8TK1KkCBqNBgBzc3MMDQ15+PChNl6sWDHu37+f8xkKIUQuS4pNZP/EvZz/Wzl0nqmlKa0me+HiVynP8xKvpvOZmIuLCxcvXgTAyMiI6tWrs3LlSpKSkoiPj2fVqlWUKVMm1xIVoqDSaDQkxiTqOw2RicfnHrGi7bIMC1hpNwd6BvWRApaP6Xwm5ufnx5w5c4iPj6dQoUKMGDGC3r17U7ZsWVQqFTExMcydOzc3cxWiwIkNi2WXfxDJ8cm8v6QzKgOVvlMS/5PWvHFo+gHFvS8DIwPcvmpM3U/qY2Aoz4LlZzoXsc8//5zPP/9c+7Ofnx9bt25lw4YNGBoa4uPjQ5MmTXIlSSEKotv7b7H9y23EPokB4NQfJ6jzcT09ZyUAoh9Gsf3Lbdw9qByJw6qcNd6z2mBfs6QeMhPZ9VoPO6dxc3OTSTCFeElyQjKHpgRz6o8T6ZYfmhqMo7sTtlVloGx9uhYYwi7/IOLD4xWxql2r4zGuBSZFTPSQmXgdb1TEhBDphV19SuDQLTy9+EQRK+5qi7G5fDjqS1JsIvsn7OX8igyaN4oWSm3eaFNRD5mJN5HtETsWL17MrVu3eP78uSKuUqkICwvLseSEKCg0Gg1nl54m+Id9pCQkp4upDFTUH9KQBkMbYWhsqKcM322Pzj0k8IsthN9Qfm45NHLEa2YbLEpa6CEz8aZ0LmITJkzg559/pmrVqnTp0gUrK6tcTEuIgiM2LJado4K4ufO6ImZR2gLvn/0o3cBBD5kJjVrDifnHOZxJ80ajkU2o83E9ad4owHQuYsuWLaNNmzYsW7Ysx3Z+8OBBfvnlF86cOcODBw/47bff6Nmzpzau0WiYPHkyixcvJjw8nLp16zJ9+nSqVKmSYzkI8SZu77vJ9q+2EftEOdNvxfcq0/IHT0yLFtJDZiLqQRTbh28l9PBdRcyqvDU+s/2wq26vh8xETtL560dMTAyenp45uvOYmBhcXV2ZPHkyZmZmivisWbP47bffmDJlCrt378bW1pYOHToQFRWVo3kIkV3J8cnsm7CH9R+uVRQwE3MTvGb64jPbTwqYnoRsvcpy78UZFrBq3WvQY0tvKWBvCZ3PxNzc3Lhw4UKO7tzLywsvLy8APv3003QxjUbDnDlzGDZsGO3btwdSx290cXFhzZo19O3bN0dzEUJXYVefEvjFFp5eUjZv2Ncuic9sP4o6WeV9YoLEmET2j9/DhVXnFLFCVqnNGxV8pXnjbaLzmdi0adMICgpi2bJl2uGnctPt27d59OgRLVu21C4zMzPD3d2do0eP5vr+hXiZRqPhzF8nWdF2maKAqQxUNBjaiC5ruksB05NHZx+ywm9phgXMwd2JnkF9pIC9hTI9E2vYsKFiWWJiIl988QWjRo2iVKlSGBqm77RSqVQcOXIkRxJ79OgRALa26efrsbW15cGDBzmyDyF0Ffs0hh0jg7i1WzmXnoWDJd4/t6F0fWne0Ad1ipoT845zZMZB1MkvNW8YG+A+ogl1BtaX0VLeUpkWseLFi6NSpf+l29raUqFChVxP6kUv56DRaBTLXhQSEpLbKenF2/q+cktOHq8nxx5xZtopEp8nKGIlW5Sm2tCaxJrHFejfUUHNPe5xHGemnOTZmaeKWBFHc2qNrYulixXXrl/L0f0W1OOlL29yvFxcXLKMZ1rEtmzZ8to7zQl2dnYAPH78GAeH//+G+/TpU8XZ2Yte9YYLopCQkLfyfeWWnDpeyfHJHJy8n9N/nlTETMxNaDHRk8odXd94P/pWUP++QrZc4dCY/SREKEfeqNazJh7fNMO4cM4/XF5Qj5e+5PbxyrcjdpQpUwY7Ozv27NlDnTp1AIiPj+fw4cNMmDBBz9mJt93TK08I/HwLYVeU3/BL1imF96w2cu9LTxJjEtk3bjcXV59XxApZm+E5xQtnbyky7wqdi9iSJUvYsWMHS5cuzTD+4Ycf4uPjQ48ePXTeeXR0NDdupN5jUKvVhIaGcvbsWaytrXF0dGTw4MHMmDEDFxcXKlSowPTp0ylSpAidO3fWeR9CZEdq88YpDgTsIyUhJV0srXmjwWduGBjJw7H68PD0AwKHbiHiVrgi5tikDF4/+WJuZ573iQm90bmILVq0iHr1Mh+B297enoULF2ariJ06dYp27dppfw4ICCAgIIDu3bszZ84chg4dSlxcHCNHjtQ+7Lxu3TosLGR4GJHzYp7EsGNEILf33lTELB0s8Z7lR6l6Mi29PqhT1Pw35xhHZx7KsHmj8aim1B5QT5o33kE6F7Hr16/Tp0+fTONVqlRh5cqV2dp506ZNCQ8PzzSuUqkYM2YMY8aMydZ2hcium7uus2NkIHFhcYpYpfer0GKiJ6aWpnrI7N2UlJxIQlIcRQpZEv0giqBhW7l3NFSxnrVzMXxm+1Gimp0eshT5gc5F7FWD+z579gy1Wp1pXIj8KDk+iQOT9nNm8SlFzMTif80bHQp+80ZBodao2XxkKdfvnycxOYEil4qjWm1JcnSyYt3qvWrS9JvmGJsZ6yFTkV/oXMRq1qzJ6tWr+eyzzyhUKP1QOnFxcaxevZoaNWrkeIJC5Janl58Q+Plmwq4qv5yVrFsK75+leSOv7TixmhNX96GJ11Bka0lSThcG0hewQtZmeE71xtkrbx/3EfmTzkXsyy+/pFOnTnh7e/Pll1/i6uqKSqXiwoULzJw5k5CQEFatWpWbuQqRIzQaDaf/PMnByfuVzRuGKhoObUT9IdK8oQ/X71/A4K4JFmsdMXymvHzr1LQMXjN8KSLNG+J/dC5iLVq04Pfff2fUqFHpxi3UaDRYWFjwyy+/5PgAwULktJjHMewYsY3b+24pYpaORfGZ5UfJuqXyPjGBOkVNcqARRbc5o1K/1KBhBE3HNKd2v7rSvCHSydZzYt26dcPPz4/du3dz69YtNBoN5cqVo2XLltIxKPK9G7uuszOT5o3KHV1pPqEVphbSvKEPkaERBA3fhuZYYV4uUSm2CTSf3oo6zTPvjhbvLp2KWFxcHI0aNWLQoEEMGjRIO6q8EAVBcnwSwT/u4+yS04qYiaUpLX/wpFJ7maNOX65svMzur3eQGKkc1ivZLRLnQeWo09xND5mJgkCnImZmZkZkZCQmJjk/hIsQuenJxccEfrGFZyHK5o1S9Uvj/XMbLB2K6iEzkRidyN7vdnJp7UVFzLioMU5D7XHr3IziRWXeL5E5nS8nenl5sX37dvr165eb+QiRIzRqDacWneDQlGBSEpXNG27D3Kk3pKFMS68nD07eJ2joFiLuRChiTh5l8ZruI80bQic6F7Hhw4fTt29fPvroI/r27Uu5cuUynI05q8F5hcgL8WHxrJ+wljv7byliRZ2K4j3Lj5J1pHlDH9Qpao7/epSjsw6hSUk/L6GhiSGNR3tQq28dad4QOsvWzM4Aly5dYuPGjZmu9+zZszfPSojXcDsyidmLLmA3PxjTGOX9lSqdq9Ls+5bSvKEnkXcjCBq+lfvH7ylixVxs8PmlLbZV5Evw2yL5+XmS766lWNRzEo0aYFymCyqDnH8wXeciNmrUqCzn8RJCn87fj+XXL3dS5fBVRczE0pRWk1pTsV1lPWQmAK5suJTavBGVqIjV7FObJmM9MCokI2+8LZLubyfx+kJIiqQQkHTzCinh5yhUa1KO1xGdi5iMXyjyqycXHrO5/waqPFDeX0msYk/fP97DsrSlHjITCVEJ7P12F5f/VTZvmNmY0XqaD+VaOeshM5FbNBoNSaEbISnyxaWoIy6S8uw/jGzq5+j+8u18YkK8ikat4dQfJzg0NRizl5o3UlQq9rSojlm3WoyUAqYXD07cJ3DoFiLvKr9clGlejtbTfShiW0QPmYlcpU6CJOXvHHUCKc9O66+ITZky5ZXrqFQqRo0a9UYJCaGLmEfRbP9qG3eCbytiYdbmrO3UiHsOxelnIZeo8po6Wc3x345wdNZhZfOGqSFNxjSj5ke15fbE28rAGIzNIeGJYrlh0Zx/HlPnIjZ58uRMYyqVCo1GI0VM5IkbO66xY2QQ8c+VI2+crFWObb51STQ1poKlIf615CwsL0XejSBw2FYe/Kds3rCpmNq8UbyyNG+8zVQqFUYlvUm6sRRSYrTLDSxcMLR1z/H96VzEnj9/rlimVqu5c+cO8+bN4+jRo6xZsyZHkxPiRUlxSQT/sJdzy84oYqaWppQf25Ijtva4hsfgXLwQ4+paYlfYUA+Zvpsu/3uRPd/uzLh546PaNBkjzRvvChPH9zEwtSHp3lbiYiIpUqIqJuX7olLl/HOZqvDwcM2rV3u1vn37YmxszPz583Nic+IFISEhuLi46DsNvXp8/hGBX2zh+XXlIxylGzrgNbONtnlDjlf2vOnxSohMYM+3O7my/pIiZla8cGrzRsvyb5JiviJ/X9mT28crxxo7mjZtyvjx43Nqc0IAqc0bJxf8x6FpwaiTXpqW3sgAty8bU3dQfRl5Q0/uH79H0LAtRIZGKmJlW5TDc5o0b4jclWNFLCQkBI0mR07qhAAg+lE027/cxt0DyuaNomWt8Jnth33NknrITKiT1Rz75TDHZh9Bo1Y2bzQd24wafaR5Q+Q+nYvYwYMHM1weERFBcHAwCxYs4P3338+pvMQ77npQCDv9t2fYvOH6QTWafd8SkyIyILU+RNwJJ2joVh6cvK+I2VQujs9sP4pXkuYNkTd0LmJt27bN8FuVRqPB0NCQTp066dSGL0RWkmIT2T9xL+f/PquImRYtRKuA1rj4VcrzvESqy+v+17wRrWzeqNWvDo39PTAqJI+firyj81/bxo0bFUVMpVJhZWWFk5OTTIop3tjjc48IHJpx84aDmyNeM32xKCUt8/qQEJnAnm92cmWDsnmjsG1hWk/3pWzzcnrITLzrdC5iTZs2zc08xDtMo9Zwcv5xDk0/kGHzRqMRjakzUJo39OXe8VCChm0lKqPmjZblaT3Nm8LFpXlD6Mcri9iSJUv47bffuHXrFsWKFaNDhw58//33MkGmyBHRD6NSmzcO3lHErMpb4zPbD7vqMimiPqiT1RyddYjjvx7NoHnDiKZfN6PGh7WkeUPoVZZFbPXq1QwdOpQiRYpQtWpV7t27x9y5c1Gr1VmO4CGELq4FhrDLP4j48HhFrGq36jQb1wLjwvJlSR8i7oQT+MUWHp56oIgVr2KLz2w/bCoW10NmQqSXZRGbP38+5cuXZ9u2bZQoUYLk5GQ++eQT/vrrL7777jsKFy6cV3mKt0hSbCL7J+zl/IpMmjcme+HSpqIeMhMajYbL/15k77e7MmzeqN2/Lu6jmkrzhsg3svxLvHTpEv7+/pQoUSJ1ZSMjhg8fzrp167h16xaurq55kqR4ezw695DAL7YQfkM5jJlDI0e8ZrbBoqQ0CelDQkQ8u7/ZydWNlxWxwrZF8JrhQ5lm0rwh8pcsi1hMTAz29unvR5QqlTqt+9OnT3Mvq/9JSUkhICCAf/75h0ePHmFnZ8cHH3zA6NGjMTKSb4IFiUat4cT84xzOqHnD2AD3EU2oM7C+TEuvJ/eOhRI0bAtR96IUsXKeznhO9aawjVx5EbpRp6j5b86x1MGey+Tuvl5ZCfR50/bnn39m4cKFzJkzB1dXVy5cuMDgwYMxMTGR0fILkKgHUWwfvpXQw3cVMWne0J8TV/fx3+V9xKxPwXC3OWjS/1s3NDXC49vmVO9VU5o3hM4i70WyffhW7h0NxayYGY3meEAuDjX5yiI2a9YsVq1apf05KSkJgPHjx1OsWLF066pUKv75558cS+7YsWP4+Pjg6+sLQJkyZfD19eXEiRM5tg+Ru0K2XmXX6O0kRCibN6r1qIHHt82leUMPToYEExS0BpOVthiHKi/fFne1xXd2W4q52OghO1FQXd18md1jdpAQmQBA3LM4zk49SbU11XLtKkuWRczBwYGIiAgiItLP0uno6MiTJ0948iT9pGc5/W3Nzc2NP/74g6tXr1KxYkUuX75McHAww4cPz9H9iJyXGJPI/vF7uLDqnCJWyKoQraZ4U8FHRgLXB41Gw9ElwZitdMAgUTlVTe0B/2veMJVL9kI3idGJ7P1+F5dWX1DEIkIiiLgTjlVZ61zZd45NxZIbNBoNP/zwAz/99BOGhoYkJyczYsQIvvnmm0xfExISkocZioyEX3nO6UkniL0Xo4jZ1C5OTf86FCpupofMRFJUIudnneXBXuWklWrzJMz6qGjVoZMeMhMFVfilZ5wOOEnsfeW/9+J1bakxqg6FbAq99vZfNY1Lvv6qtW7dOlauXMnChQupXLky586dY/To0Tg5OfHhhx9m+Jq3cZ6fgjJ/kTpFzYl5xzky4yDq5AyaN0Y2pc7H9XK9eaOgHK+8du/oXQKH7Sb6vrJ5I6FSJLHv36dxqx5y7F5B/r5SqVPU/Pf7MY7MPIgm5aWH4U0McfdvSu1+dbl2/VrBmE8sN3z33Xd89tlndOqU+s2watWq3L17l5kzZ2ZaxIR+RN2PZPvwbYQeUTZvWDsXw2eWHyWq2+khM5GSlMLRmYc4/vtReOm6i8ZYTYz3AxLrR1DGviJ1XTz0k6QoUCJDIwgavo37x0IVsWIViuEzuy22VUvkSS75uojFxsZiaJj+mr2hoSFqtTqTVwh9CNlyhV1jdmTcvNGzZmrzhplMS68P4beeE/jFFh6deaiI2VSxoehgUx4ZaKhd6X1qVWiMoUG+/kgQ+cCVjZfZ/fUOEv/XvPGiGr1r0eTrZnn67z1f/8X6+Pjw888/U6ZMGSpXrszZs2f57bff6Natm75TE6Q2b+wbt5uLq88rYoWszfCc6o2zVwU9ZCY0Gg0XV59n37jdJMUmKeJ1Btaj0YgmGJkayeUxoZPE6ET2freTS2svKmJmxczwnO5D+VbOeZ5Xvi5iU6dO5ccff+Srr77i6dOn2NnZ0adPH3lGLB94eOYBgV9sIeJWuCLm1LQMrWf4Ym5nnveJCeIj4tk9ZjshW64qYkVKFMFrZhucmuTyE6jirfLg5H2Chm4h4k6EIlamWVlaT/elSAn9zGSQr4uYhYUFkydPlsGG86HoB9GKAmZoYoj7qKbU7l9X5+aN+GQNi6/GcPJpIs1LmtLFuTBGMmrHaws9fIeg4duIfqBs3ijvVQHPKV6YFZORN4Ru1Clqjv96lKOzDmXYvNF4jAe1Pqqj15F2XquIRUREEBqaekPPwcGBokWL5mhSIv+r4ONCte41tIP4vs7N3GfxKbwf9JQLz5JJAf69GceSkFjWeRXHzEgKWXakJKVw5KeD/DfnmKJ5w6iQER7jWlCtew0ZeUPoLDI0gqBhW7l/XPk4RjEXG3x+aYttFVs9ZJZetorY4cOHGT9+PMeOHUu3vEGDBnz33Xe4u7vnaHIif/P4rjmhR+/i6O5E02+y37zx7X+RnH2WrP05UQ2HHyUy+1wU/rVlBmddPb/5nMAvNvP47CNFrEQ1O7xntaFYBRl5Q+guq+aNmn1q02SsB0aF8kezls5FbPv27fTs2RNzc3P69+9PhQoV0Gg0XL9+nTVr1tC+fXuWL1+Ol5dXbuYr8hHjwiZ039QbE/PXGzYqJFzZcABw/IlyChChpNFouPjPefZ9n0HzhgrqflKfRl81wdBEOSqHEBlJiEpg73e7uLwug+YNGzNaT/OhnB6aN7KicxEbP3485cqVIygoCGvr9MOHjBkzBi8vL8aPHy9F7B3zugUMoLCRQYbLLU3kkterxIfHsWv0dq5tU45QU8TOHO+ZbXBs7KSHzERB9eDEfQKHbiHybgbNG83L0Xq6D0Vs9dO8kZWMP0UycP36dfr06aMoYADFihWjT58+XL9+PUeTE2+3gVWKUMw0fcGyNzPgy+oyn1hW7h66w3LvxRkWMGdvF3oG9ZECJnSmTlZzdNYhVndZoShghqaGNPu+Je3/6pgvCxhk40ysbNmyxMQox8ZKExMTQ5ky0rYrdNemjBkTE9UsuhxDeIKa4mYGfFnDgmo2Mqp9RlISUzj800FOzM2gecPMiGbjWlK1W3Vp3hA6i7wbQeCwrTz4T9m8YVOpOD6z/VLnBMvHdC5i/v7+fPXVV7Rs2ZJ69eqlix0/fpwFCxYwY8aMHE9QvN16uhShp0v+/IaXnzy//ozAoVt4fC6D5o3qdvjM8sPauVgGrxQiY5f/vcieb3eSGKW8B12rbx0aj26ab5o3sqJzETtw4AD29vZ4eXlRu3ZtnJ1Tb+5dv36dU6dOUaVKFYKDgwkODta+RqVSMX369JzPWoh3hEaj4cLKc+wbv5vkuOT0QRXUHdSARl82luYNobOEyAT2fLuTK+svKWJmxQvjNd2Hsi3K6yGz16NzEVu0aJH2/0+ePMnJkyfTxS9evMjFi+k7WqSIvb0SUjQ8jkvBvrAhxvJwcq6Iex7HLv/tXA9S3vsytzfHa2YbHN3l3pfQ3f3/7hE0dAuRoZGKWNkWqc0bZlagjn+CyrR4gbg0rXMRe/78eW7mIQqQiSci2HArjucJamwKGdKjQmGG1ZBmjJx09+Adtn+5leiH0YpYBV8XWk32opCVzMkmdKNOVnPsl8Mcm30EjfqlkTdMDWk6thnVe1Ul8dJ04q5cRqNORGVaHBPnfhjZ1NVT1rrJ18NOifxn6dUY5lyMJvZ/V7bCEpKZeTYKV2sjvBzlQ/VNpSSmcHj6AU7MP65o3jAubEyz71vi+kG1AvENWeQPEXfCCRq6lQcn7ytiNpX/17xRyZaEi9NIefL/t4M0SREkXvkFw/q/oDLOv19Ss13Erl+/zo4dO7h7N3XeKEdHRzw9PalQQUYrfxesvhGrLWBpIpI0LLoSK0XsDT27FkbQ0K08Pp9B80YNO3xmt8W6XO5M8S7eTpfX/a95IzqD5o1+dWjs74FRISM0Gg0pEZcV62jiH5J0bwsmZfPvzCE6F7Hk5GRGjhzJkiVLFPN5jR07lt69ezNjxgyMjOTk7m2WnMlUbkkvXaIQutNoNJxfcZb94/eQHK9s3qj3aUPchrtjaCzNG0I3CZEJ7PlmJ1c2KJs3CtsWpvV0X8o2L/fCUg0aTbJiXQCSlYNJ5yc6V5yJEyfy119/0a1bNz755BNtd+K1a9eYN28eS5YswcrKivHjx+daskL/GpQw4fCjxHRXuoxU0Lp0Ib3lVJDFPYtl1+jtXA+6poiZl7TAe6YvDo2keUPo7t7xUIKGbSUqg+aNcq3K4znVm8LF0z/WolIZYGBWGnX8S1cBjK0wLNUmN9N9YzoXsRUrVvD+++8zZ86cdMtr167N3LlziYuL4++//5Yi9pYbW9uS88+SOPY4kcgkDdYmKpqWNOXjKvKsV3bdOXCb7cO3EvNYOYhAhTYVaRXQWpo3hM7SRt44/uvRDJo3jGj6TTNq9K6V6f1Uk8rDSDw3HnXMHdAkoTIpjlHpNhgWLp0X6b82nYtYbGwsTZo0yTTu4eHB7t27cyQpkX+ZGKpY41Wck08SOfE0kcb2Jrhaywgb2ZGckMzhaQc4ueA/Rcy4sDHNxrfEtYs0bwjdhd8OJ2joFh6eeqCIFa9ii89sP2wqFs9yG4ZmJShUfzYpYf+hSQjDyLYRKpP8fw9W5yLm7u7OkSNH6N+/f4bxI0eOyFQs75A6tibUsZXilV3ProUR+PkWnlx8rIjZ1bTHZ7YfVmXz/weHyB80Gg2X111k73e7MmzeqN2/Lu7+TTEy1e2jXqUyxKh4w5xOM1fpXMRmzJhB586d+eqrrxg4cCDly6c+0X3jxg3mzZvHuXPnWLNmTa4lKkRBptFoOP/3WfZPyLh5o/4QNxoOayTNG0JnCRHx7P5mJ1c3KrsKC9sWwWuGD2WalcvglW8XnYtY/fr10Wg0XL16lT///FN7qUOjSb32amRkRP369dO9RqVScf++8tkEId4lcc9i2em/nRvbM2jeKGWBz89tKN3QUQ+ZiYLq3rFQgoZtIeqesnOwXKuyeE5rQ2GbwnrILO/pXMQ6dOgg1+iFyKbbwbfYPnwbsU+UzRsubSvRclJrChWVzk6hm5SkFI7NOszx3zJo3jBOoWG3q1R9/wFmxTroKcO8p3MRe7krUQiRueSEZA5NDebUwhOKmHERY5pPaEWVTlXli6HQWVbNG8Uco2gx6BzWpWLRRD0m5VEwRvbN8z5JPZAnk4XIYWFXnxI4dAtPLz5RxOxq2eMzS5o3hO40Gg2X1lxg77hdJMUkKeLVvW9Tr9M1DI3/d2amSSIl/KwUscw8ePCAM2fOEBERoRi5A6B79+45kpgQBY1Go+Hs0tME/7CPlIT0zRsqAxX1hzSkwVBp3hC6i4+IZ/fYHYRsvqKIFSlRBI9PblC6/EuzHKiMMbCqlkcZ6p/ORSwxMZHPPvuMtWvXolarUalU2qaOFy+JSBET76LYsFh2jgri5s7riphFaQu8f/ajdAMHPWQmCqrQI3cJGr6V6PvK5o3yrZ3xnOqNUfwBEq/dhuT/H53DwLICRiWa5WWqeqVzEZs0aRJr165lzJgxNGrUiLZt2zJnzhzs7e359ddfefLkCXPnzs3NXIXIl27vv8X2LzNu3qj4XmVa/uCJqTRvCB2lJKVwdOYhjv9+VDGTgVEhIzy+a0G1HjX+d/LghcrUhqTQ9ZCSgIFFeUzK9UFl8O6c7etcxNauXUvXrl0ZMWIEz549A6BkyZI0a9aMZs2a0aZNGxYtWsS0adNyLVkh8pPk+GQOTg3m9B/K5g0TcxOaT2xF5Q6u0rwhdBZ+6zmBX2zh0ZmHipitawl8ZvtRzMUm3XIjm7r5fs6v3GSg64qPHz+mYcPUJ7nTRqqPj48HUi8ntm/fno0bN+ZCikLkP2FXn7Lq/eUZFjD72iXpse1DqnSU7kOhG41Gw4V/zvG375IMC1idgfX4YH0PRQET2TgTs7GxITw8HAALCwvMzMy4deuWNp6UlERMjPJyypt6+PAh33//PTt27CA6OpqyZcsyY8aMLMdxFCK3aDQazi45TfCPmTRvfNaQBl9I84bQXXxEPLvHbCdky1VFrIidOV4/+eLUpIweMisYdC5i1atX5/jx40DqmVfjxo35/fffqVGjBmq1mvnz51O9evUcTS48PBxvb2/c3Nz4559/sLGx4fbt29ja2ubofoTQRezTGHaMDOLW7huKmIWDJd4/t6F0fWneELoLPXyHoOHbiH6gbN5w9q5AqynemFnLTAZZ0bmIffTRRyxbtoz4+HgKFSrExIkTadu2LX5+fmg0GooVK8aPP/6Yo8nNnj0be3t75s2bp11WtmzZHN2HELq4tfcmO0ZsI/ZJrCImzRsiu1ISUzgy8yD/zTmWcfPGuBZU615DLkfrQOci5uvri6+vr/bnypUrc/LkSYKDgzE0NMTNzQ0rK6scTW7Lli20atWKvn37EhwcjL29PR9++CEff/yx/HJFnkiOT+bg5P2c/vOkImZibkKLiZ5U7uiqh8xEQfX8xjMCh27h8dlHiliJanZ4z2pDsQpy70tXqvDw8Hw7r7ydnR0An376Ke+//z7nzp3D39+fcePGMXDgwAxfExISkuFyIbIr6mYkpyedIOqmcoZcK1drao2pS+GSMhmo0I1GoyE08A4XfztHSnxK+qAKyn9QgYofVcHAWOd+u3eCi4tLlvHXLmLh4eH07t2bH374gZo1a75Wcq9ia2tL7dq12b59u3bZhAkT2Lx5M8eOHcuVfeZHISEhr/xFiv/3psdLo9FwdvEpgiftIyUh/YeNykBFg6GNaPCZGwZGb8eHjfx9Zc/rHK/48Dh2jd7OtW3KL9lF7MzxntkGx8ZOOZVivpLbf1+vPXZiYmIiBw4c0HYs5gY7OzsqVaqUblnFihUJDQ3NtX2Kd1vMkxh2jgzk1p6bipilgyXes/woVS9/T9cu8pe7h+6wffhWoh9GK2LO3i60muIlzRtvIF8PAOzm5sa1a+nnYLp27RqOjjL3ksh5N3ffYMfIQOKeKps3Kr1fhRYTPTG1NNVDZqIgSklM4fBPBzkxN4PmDTMjmo1rSdVu1eX+/hvK10Xs008/xcvLi+nTp9OxY0fOnj3L/Pnz+fbbb/WdmniLJMcncSBgP2f+OqWImVj8r3mjgzRvCN09v/6/5o1zGTRvVLfDZ5Yf1s7F9JDZ2+e1i1ihQoXo3r07JUuWzMl80qlTpw7Lly9nwoQJTJs2DQcHB8aOHcuAAQNybZ/i3fL08hMCv9hC2JWniljJuqXw/rkNRZ2s8j4xUSBpNBourDzHvvG7SY5L/zA8Kqg7qAGNvmyMoYk8DJ9TXruIWVpa8vvvv+dkLhny9vbG29s71/cj3i0ajYbTf57k4OT9yuYNQxUNhzai/pC3p3lD5L6453Hs8t/O9SBl84a5vTleM9vg6P52Nm/o02sVsZiYGJ4/f66diuVFcr9K5Hcxj2PYMWIbt/fdUsQsHYviM8uPknVL5X1iosC6c+A227/cRswjZfNGBV8XWk32opCVNG/khmzNJzZ16lQWL15MWFhYpuuljXAvRH50c9f11OaNsDhFrHJHV5pPaIWphTRvCN2kJKZwePoBTsw/rmjeMC5sTLPvW+L6QTVp3shFOhexUaNGsWTJEnx9fWncuHGOj84hRG5Kjk8i+Md9nF1yWhEzsTCh5Y+tqdS+St4nJgqsZ9fCCPxiC08uPFbEStSww2d2W6zLWeshs3eLzkVs/fr19OjRg19//TU38xEixz259ITAzzfzLER5BaFU/dJ4z2yDpWNRPWQmCiKNRsO5v8+wf/wekuOVzRv1BjfAbbg0b+QVnYuYWq2mXr16uZmLEDlKo36heSMxo+YNd+oPaSjNG0Jncc9iOfn9MR4dVM75ZV7SAu+Zvjg0kuaNvKRzEWvVqhVHjhzho48+ysV0hMgZMY+i2TEyMMPmjaJORfGe5UfJOtK8IbJn77jdGRawCm0q0iqgtTRv6IHOX0GnTp3K+fPnmTRpEo8eKR/gEyK/eHT4Ict9FmdYwKp0rkr3rR9KAROvpcnYZhhbGGt/Ni5sjOc0b9r83k4KmJ7oPACwvb09Go2GpKQkAIyNjTEwSF8DVSoV9+/fz/ks33EyQKvuTi78j+CJexXLTSxNaTWpNRXbVc7znPI7+fvKnuBF+zg5/jh2Ne3xme2HVVlp3shKvhkAuEOHDtImKvK1xKQEwhzugIkGEv//b7VUAwe8Z/pi6SDNG+L1aRIjSLqzlspVrmM3tSHlO/hiZGL86heKXKVzEZszZ05u5iHEG4mOi+Cv7VN59DwU0zbWWKx3AAMNDYa50/CzRhgYSvOGeH0pUddIOPcjmvgHFAYK250j+fIpDKt/J1/u9SxfDwAshK62HfubR89Tp+hJqPUcw8emJFaNIKF5RSlg4o0lhixAE//g/xeoE0l5dpKUsOMYFW+gv8RE1kXsxIkT2d5g3bp1XzsZIV7Xs+gn//+DCmK9UzvI7j6+rqeMxNtEk6AcIBp1AsmP90sR07Msi5inp6fOp8oajQaVSiXDTgm9MDHKeKgoM9MieZyJeCsZZtx5aGAmXa76lmUR++233/IqDyHeiFvl1jx4doe4hP8fgNWysDXNarbTY1bibWFUoilJsaGgjtcuU5k5YOzYXo9ZCXhFEevRo0de5SHEG6lSpg4JyfEcvbyTyOhwilkWp2WtDtgWlW/K4s0Zl/kAUJH85ACJcRGYWjpg4vIJKiM509c3aewQb41azu7UcnaX555EjlOpVJiU/QCTsh9wT/6+8hVp2xJCCFFgSRETQghRYEkRE0IIUWBJERNCCFFgSRETQghRYOk8ir0QQgiR38iZmBBCiAJLipgQQogCS4qYEEKIAkuKmBBCiAJLipgQQogCS4pYPvXTTz/RokULHB0dcXZ2pmvXrly8eFHfaRUYM2bMwMrKipEjR+o7lXzr4cOHDBo0CGdnZ+zs7GjYsCEHDhzQd1r5UkpKCj/88AM1atTAzs6OGjVq8MMPP5CcnKzv1PKFgwcP0q1bN6pUqYKVlRXLly9PF9doNAQEBFC5cmXs7e3x8/Pj0qVLObJvKWL51IEDB+jfvz9BQUFs3LgRIyMj3n//fZ4/f67v1PK948ePs3jxYqpWrarvVPKt8PBwvL290Wg0/PPPPxw9epSpU6dia2ur79TypZ9//pmFCxcyZcoUjh07xuTJk1mwYAE//fSTvlPLF2JiYnB1dWXy5MmYmSnnXps1axa//fYbU6ZMYffu3dja2tKhQweioqLeeN/ynFgBER0djZOTE8uXL8fX11ff6eRbERERNGvWjFmzZjF16lRcXV2ZNm2avtPKdyZMmMDBgwcJCgrSdyoFQteuXbG2tmbu3LnaZYMGDeL58+esWrVKj5nlP6VLl2bq1Kn07NkTSD0Lq1y5Mh9//DEjRowAIC4uDhcXFyZOnEjfvn3faH9yJlZAREdHo1arsbKy0ncq+dqwYcNo3749zZo103cq+dqWLVuoW7cuffv2pUKFCjRp0oT58+ej0ch32oy4ublx4MABrl69CsDly5cJDg6mdevWes4s/7t9+zaPHj2iZcuW2mVmZma4u7tz9OjRN96+zCdWQIwePZrq1avToEEDfaeSby1evJgbN24wb948faeS7926dYs//viDTz/9lGHDhnHu3Dn8/f0BGDhwoJ6zy3+GDRtGdHQ0DRs2xNDQkOTkZEaMGMGAAQP0nVq+9+jRIwDFpWpbW1sePHjwxtuXIlYAjB07liNHjhAYGIihoaG+08mXQkJCmDBhAtu2bcPExETf6eR7arWa2rVrM27cOABq1qzJjRs3WLhwoRSxDKxbt46VK1eycOFCKleuzLlz5xg9ejROTk58+OGH+k6vQFCpVOl+1mg0imWvQ4pYPjdmzBjWrVvHpk2bKFu2rL7TybeOHTtGWFgYjRo10i5LSUnh0KFDLFq0iPv372NqaqrHDPMXOzs7KlWqlG5ZxYoVCQ0N1VNG+dt3333HZ599RqdOnQCoWrUqd+/eZebMmVLEXsHOzg6Ax48f4+DgoF3+9OnTHGkkkiKWj/n7+7Nu3To2b95MxYoV9Z1Ovubn50ft2rXTLRsyZAjOzs58+eWXcnb2Ejc3N65du5Zu2bVr13B0dNRTRvlbbGys4iqIoaEharVaTxkVHGXKlMHOzo49e/ZQp04dAOLj4zl8+DATJkx44+1LEcunRowYwapVq1i2bBlWVlba68pFihTB3Nxcz9nlP1ZWVoqml8KFC2NtbY2rq6t+ksrHPv30U7y8vJg+fTodO3bk7NmzzJ8/n2+//VbfqeVLPj4+/Pzzz5QpU4bKlStz9uxZfvvtN7p166bv1PKF6Ohobty4AaReqg4NDeXs2bNYW1vj6OjI4MGDmTFjBi4uLlSoUIHp06dTpEgROnfu/Mb7lhb7fCqzLkR/f3/GjBmTt8kUUH5+ftJin4WgoCAmTJjAtWvXcHBw4OOPP+aTTz7JkfsUb5uoqCh+/PFHNm/ezNOnT7Gzs6NTp06MGjWKQoUK6Ts9vQsODqZdu3aK5d27d2fOnDloNBomT57MX3/9RXh4OHXr1mX69Ok58gVTipgQQogCS54TE0IIUWBJERNCCFFgSRETQghRYEkRE0IIUWBJERNCCFFgSRETQghRYEkRE/mGn58f9evXz/X9BAQEyGwAOWz06NF4eXmlW2ZlZUVAQICeMtKvwMBAHBwcePr0qb5TeetJERNay5cv1458sX///gzXadmyJVZWVnlSbN7EvXv3CAgI4OzZs/pOJUMJCQnMmzcPLy8vnJycsLW1pUaNGgwZMoTTp0/nyj4PHz5MQEAA4eHhObrd0NBQ/vzzT+1cUblpxowZ9OjRg8qVK2NlZcXw4cMzXTcqKoqRI0dSsWJF7O3t8fT0ZPfu3Rmue+PGDbp3746TkxOOjo50796dmzdvZrju7t278fT0xN7enooVKzJq1Ciio6PTrePt7Y2Tk5NMmpkHpIgJhUKFCrF69WrF8uvXr3Py5MkCMULB/fv3mTJlCufOnVPERo4cycOHD/WQVaqwsDB8fX3x9/fH0tISf39/ZsyYQZcuXThy5AgtWrTg3r17Ob7fI0eOMGXKFCIiInJ0u/PmzaNEiRKKubUePnzIyJEjc3RfEydO5Pjx49SsWTPL9TQaDT179mTp0qX07t2bgIAADAwM6NKlC/v27Uu37qNHj/D19eXcuXOMGjUKf39/zp07h6+vL48fP0637r59++jSpQuGhoYEBATQu3dvlixZQo8ePdLNxaZSqfjoo4/466+/iIyMzLkDIBRk7ESh4OXlxYYNG5g+fXq6kd9XrVpFiRIlcHZ2ztHLJLGxsRQuXDjHtvcqRkZGGBnp709/8ODBnDp1ikWLFtGxY8d0sbFjx/Lrr7/qKbPsS05OZtWqVXTv3l0xXFVufNk5ffq0djaHrC4Jb968mf379zN37lzt+IY9evSgcePGfP311xw4cEC77k8//cTz5885evQo5cqVA8DX1xc3Nzd++uknJk+erF137NixVKhQgU2bNmkHlS5fvjxDhgxhy5YttG3bVrtu+/btGT16NOvXr5eR7nORnIkJhU6dOhEdHU1gYGC65WvWrKFjx44YGCj/bJYvX0779u2pWLEiJUqUoG7duvz888+KUb7T7nudP3+edu3aUapUKb766qtMczl8+DBOTk506dKF+Ph4IPWMsF+/fjg7O1OiRAnc3d1ZtmyZ9jUvzrg7ZMgQ7SXStPszGd0Tq169Op06deLEiRP4+Phgb29P1apV+f333xU5hYaG0qtXL0qXLk25cuX4/PPPOX/+PFZWVixfvjyLIwv//fcf27dvp1evXooCBqkjow8dOpTSpUtrl128eJFu3brh5OREyZIlad26NTt27FC8duHChbi7u1OqVCnKli1Ls2bNWLRokfY9jx8/HkidOyztmAQHBwOpxaFLly44Oztjb29PzZo1+eSTT4iJicny/Rw9epTHjx/TvHlzRezle2Jpl6sPHTrEhAkTqFSpEvb29nTo0IFbt25luZ80uk5H9O+//2JtbU2XLl20y0xNTenTpw/nz58nJCREu3z9+vV4enpqCxiAs7MzLVu25N9//9Uuu3r1KhcuXKBPnz7pZkXo2rUrRYsWTbcupE5BUqVKFTZv3qxTzuL1yJmYUChVqhSNGzdm9erVtG/fHkj98L1x4wYffPBBhpfoFixYgIuLC56enpiZmbFnzx6+//57IiMj+e6779KtGxERQceOHWnXrh2dOnWiaNGiGeaxb98+unfvTsuWLVm0aBEmJiZcuXIFb29vbGxsGDJkCEWLFmX79u189tlnREZG8umnn1KpUiVGjx7N5MmT+eijj7RzjFWtWjXL93379m26detGjx496NKlC+vWrWPs2LFUrlxZO7V6bGws7733HqGhoQwcOBAnJyc2b97M4MGDdTq227ZtA9B59PNr167h4+ODiYkJn376KUWKFOHvv/+ma9euLF68WDvo6pIlSxgxYgTvvfceH3/8MUlJSVy+fJkjR47Qr18/2rVrR0hICOvWrWPSpEnY2NgAUKlSJZ4+fUqHDh2wsbFh6NChWFlZERoayrZt24iJiaFIkSKZ5nfkyBEAatWqpdP7gdSzGTMzM4YPH05YWBi//vorAwcOZPv27Tpv41XOnDlDrVq1FNOn1KtXTxt3cXHhwYMHPHr0iLp16yq2Ua9ePYKCgnj48CH29vacOXMGQLGukZERtWrV0sZfVKdOHTZu3JhjE0AKJSliIkNdunRhxIgRhIeHY2VlxapVq3B2dtbOB/SyrVu3prskOGDAAD7//HPmzZuHv79/usuSjx8/ZvLkyQwaNCjT/QcFBdGnTx/atm3L3LlztZf/Ro8erZ2bKG1//fv3p2/fvgQEBNCnTx9KlChBq1atmDx5MvXr16dr1646vedr166xfv167VlFr169qFatGosXL9YWsT///JMbN26kuxTYv39/bbF/lStXrgCvLqhpJkyYQGxsLDt37tTOKdenTx/c3d0ZM2YMfn5+GBgYEBQURJUqVViyZEmG26lWrRrVq1dn3bp1+Pn5UaZMGW1sy5YtPH/+nHXr1qWbk23s2LGvzO/q1atYWlpibW2t0/uB1ClyNm/erD2jt7a2ZuzYsVy6dIkqVarovJ2sPHz4MMPmo5IlSwLw4MED7Xrw/xM3vsje3l67jr29/SvXPXHihGJ52bJliYiI4OHDh9p9i5wllxNFhtq3b49KpWLDhg0kJyezfv36dJdmXpZWUFJSUggPDycsLIwmTZoQExOT7tINpH5z/eijjzLd1oYNG+jVqxedO3dm/vz52gIWHh7O3r17ef/994mLiyMsLEz7n6enJ1FRUZw6deq137Ozs3O6y2KmpqbUq1cv3aWunTt3UqJECd5//33tMkNDQz7++GOd9hEVFQWAhYXFK9dNSUlh165d+Pj4pJsU1dLSkn79+hEaGsqFCxe027t3716GH6SvkpZLYGAgSUlJ2Xrts2fPsv24Qt++fdNdkm7cuDGAzpcUdREXF5fhTN5py9IuTcfFxaVb/qK0e3pp67xq3bT4i9KKe1hYWLbfg9CNFDGRoaJFi+Ll5cU///zDnj17ePLkSZZF7PDhw/j6+lKyZEnKli2Ls7Mzn3zyCYCiG87e3j7Tm/6hoaH069cPLy8vfvnll3QfdtevX0ej0TBlyhScnZ3T/TdkyBCAN2o4yWhWYysrK54/f679+e7du5QrV05xX9DZ2VmnfaQVjLRilpWnT58SExOT4azelSpVAuDOnTsADBs2DHNzc1q1akWtWrUYPny4ogsvM02bNqVdu3ZMmTKF8uXL07VrV/766y9F23hmXuzK08XLxzmtCL54nN+UmZkZCQkJiuVpy9L+/szMzNItf1FaoUtb51XrpsVflHZs5FJi7pHLiSJTXbp0oU+fPkDqfYDMPqhv3bpFhw4dKF++PAEBATg4OGBqasqZM2cYN26corkjo3/saWxtbXFycmLPnj0cOXJEez8L0G4nbVbijLzJJHsv3z9Jo8uHtK4f5JUqVWLz5s1cvHgRd3f3bOWX1f4qV67M8ePH2blzJ7t27SIoKIg///yTvn37MnPmzCy3pVKpWLp0KSdOnCAwMJC9e/cybNgwZsyYwa5duyhRokSmry1WrFi2W/bf5Djrys7OTjsb+ovSLiOmXdpLuzSY0bpplw/TLiu+uO7LDSZplxxflvZMXto9SJHz5ExMZMrb2xtLS0sOHjyY5VnY1q1biY+PZ+XKlfTv3x9vb2+aN2/+WqNimJqasnLlSqpWrUrXrl3TXR5M++AwMjKiefPmGf6X9oGbW998HR0duXnzpqIwp03N/iq+vr4ArFy58pXrFi9enCJFinD16lVFLO0SrZOTk3ZZkSJFaN++PbNnz+bs2bN06dKFP//8k/v37wOvPiZ169bl66+/ZseOHaxevZq7d+9meo8tTaVKlYiMjOTZs2evfD95Ka3RIiUlJd3y//77D0D7nFmpUqUoUaJEhpdh//vvP+zs7LTFKa155eV1k5OTOXPmTIbPrt28eZOiRYtmeB9N5AwpYiJTpqamzJgxA39/fzp37pzpemnfrF/8Jp2QkMD8+fNfa7/m5uasXr2aMmXK0LFjRy5evAiknqV5eHjw119/ERoaqnjdi5cS0+7R5fToFJ6enjx+/Jj169drl6WkpLBgwQKdXl+vXj28vLxYtmwZGzZsUMRTUlKYPXs29+7dw9DQkFatWhEUFMS1a9e060RFRfHnn3/i4OCgbRB5uYgYGRlpY2nHILNjEh4erjgLSvtAftXxa9iwIUCujTLyutq3b8+zZ89Ys2aNdllCQgKLFy/G1dUVFxeXdOvu3Lkz3Qgd169fZ/fu3ekadipWrIirqyuLFy8mMTFRu3zVqlWEh4dn2Nxz8uRJGjRoIJcTc5FcThRZyqp4pWnVqhUmJiZ069aNjz76iMTERFauXJnh82S6srKy4t9//8XPz48OHTqwdetWnJ2d+emnn/D29qZx48b06dMHZ2dnwsLCOHPmDLt37+bu3btA6j0qS0tLFi1ahLm5Oebm5lSpUuWNLjcCfPTRRyxYsIDBgwdz8uRJbYt92qgMunxYzZkzh86dO9OnTx+8vLxo3rw5FhYW3L59mw0bNnD9+nXtcf/222/Zu3cvvr6+DBgwQNtiHxoayl9//aU9xh06dMDW1hY3NzdKlCjBzZs3mT9/Pq6urlSuXBlA23k4ceJEOnXqhImJCR4eHqxevZqFCxfStm1bypUrR1xcHMuXL8fQ0PCVXZcNGjTA1taWPXv2aDs4c9PKlSu1v2NIbZWfNm0akPq8VtqZ6XvvvUfjxo0ZOnQoISEhODg4sGLFCm7evMnatWvTbfOrr75iw4YNtGvXTvuoxJw5c7C2tubLL79Mt+6PP/5Ip06deO+99+jWrRuhoaH8+uuvNGnSRPu4Q5qHDx9y+fJlBg4cmOPHQfw/KWLijVWoUIHly5czYcIExo0bh42NDd26daNJkyZ06NDhtbdbvHhx1q9fT5s2bWjfvj1bt26lQoUK7N27l6lTp7J69WqePn2KjY0NlSpVYuLEidrXmpqaMm/ePCZOnMiIESNISkrC39//jYtYkSJF2LRpE/7+/tpn19q1a8fXX3+Nt7e3TqNU2NjYEBgYyJ9//snatWuZPHkycXFxlCxZUvuAcqlSpQBwcXEhMDCQ8ePH89tvv5GYmEj16tVZuXJluvuCffv2ZfXq1cyZM4eoqCjs7e3p2bMnI0eO1Ba6+vXr88033/DXX38xZMgQ1Go1mzZtonHjxpw6dYp///2Xx48fY2FhQY0aNZg6deorx8g0Njama9eu/Pvvv0yYMCHXzziWLl3KwYMHtT+fPHmSkydPAuDm5qYtYiqVihUrVjBx4kQWL15MVFQUrq6urFq1SvFgtr29PVu3buXrr7/Wjs7h7u7OpEmTFPe5WrRowT///MOkSZMYPXo0FhYW9O7dm++++07x3jdu3Iipqekb/RsQr6YKDw/PubupQryjNm3aRO/evQkMDMTNzU3f6eSpu3fvUq9ePRYvXoyPj4++08kXNBoNjRs3xsPDI92wVSLnyT0xIbLp5eeBUlJSmDt3LpaWltkaueJt4ejoSN++fWXE9hcEBQVx+/btLIdUEzlDzsSEyKZOnTpRokQJateuTXx8PBs2bODEiROMHz+eoUOH6js9Id4pUsSEyKY5c+awdOlS7ty5Q1JSEs7Oznz88cf07dtX36kJ8c6RIiaEEKLAkntiQgghCiwpYkIIIQosKWJCCCEKLCliQgghCiwpYkIIIQosKWJCCCEKrP8DcVxoBt4pDIYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fe = smf.ols(\"purchase ~ mkt_costs + C(city)\", data=toy_panel).fit()\n",
    "\n",
    "fe_toy = toy_panel.assign(y_hat = fe.fittedvalues)\n",
    "\n",
    "plt.scatter(toy_panel.mkt_costs, toy_panel.purchase, c=toy_panel.city)\n",
    "for city in fe_toy[\"city\"].unique():\n",
    "    plot_df = fe_toy.query(f\"city=='{city}'\")\n",
    "    plt.plot(plot_df.mkt_costs, plot_df.y_hat, c=\"C5\")\n",
    "\n",
    "plt.title(\"Fixed Effect Model\")\n",
    "plt.xlabel(\"Marketing Costs (in 1000)\")\n",
    "plt.ylabel(\"In-app Purchase (in 1000)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65854efd",
   "metadata": {},
   "source": [
    "花一点时间来理解上面的图片到底告诉你固定效果在做什么。请注意，固定效应模型**对每个城市拟合一条回归线**。还要注意拟合的回归线是平行的。线的斜率是营销成本对应用内购买的影响。所以**固定效应是假设因果效应在所有实体中都是常数**，在我们的例子中，实体是城市。这可能是一个弱点或一个优势，取决于你的目的。如果您的兴趣是对每个城市找到其包含的因果关系，那么这是一个弱点。由于 FE 模型假设这种效应在实体之间是一样的，因此您不会发现因果效应有任何差异。然而，如果你想找到营销对应用内购买的整体影响，面板结构的数据是一个非常有用的途径来使用固定效应模型进行探索。\n",
    "\n",
    "## 时间效应\n",
    "\n",
    "就像我们在个人颗粒度做固定效果一样，我们可以为时间颗粒度设计一个固定效果。如果为控制固定的个体特征而为每个个体添加一个虚拟变量，则添加一个时间颗粒度的虚拟变量将控制每个时间段内固定的变量，但不同时间之间是可以变化的。这种变量的一个例子是通货膨胀。价格和工资往往会随着时间的推移而上涨，但每个时间段的通货膨胀对于所有实体来说都是相同的。举一个更具体的例子，假设婚姻随着时间的推移而增加。如果工资和婚姻比例也随时间变化，我们就有时间作为一个混淆因素。由于通货膨胀也使工资随着时间的推移而增加，我们在婚姻和工资之间看到的一些正相关可能仅仅是因为两者都随着时间的推移而增加。为了纠正这一点，我们可以为每个时间段添加一个虚拟变量。在“线性模型”中，这就像在我们的公式中添加“TimeEffects”并将“cluster_time”设置为 true 一样简单。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5c000d3d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Parameter Estimates</caption>\n",
       "<tr>\n",
       "     <td></td>     <th>Parameter</th> <th>Std. Err.</th> <th>T-stat</th>  <th>P-value</th> <th>Lower CI</th>  <th>Upper CI</th> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>expersq</th>  <td>-0.0062</td>   <td>0.0008</td>   <td>-8.1479</td> <td>0.0000</td>   <td>-0.0077</td>   <td>-0.0047</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>hours</th>    <td>-0.0001</td>  <td>3.546e-05</td> <td>-3.8258</td> <td>0.0001</td>   <td>-0.0002</td> <td>-6.614e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>married</th>  <td>0.0476</td>    <td>0.0177</td>   <td>2.6906</td>  <td>0.0072</td>   <td>0.0129</td>    <td>0.0823</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>union</th>    <td>0.0727</td>    <td>0.0228</td>   <td>3.1858</td>  <td>0.0015</td>   <td>0.0279</td>    <td>0.1174</td>  \n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mod = PanelOLS.from_formula(\"lwage ~ expersq+union+married+hours+EntityEffects+TimeEffects\",\n",
    "                            data=data.set_index([\"nr\", \"year\"]))\n",
    "\n",
    "result = mod.fit(cov_type='clustered', cluster_entity=True, cluster_time=True)\n",
    "result.summary.tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e0db7ec",
   "metadata": {},
   "source": [
    "在这个新模型中，婚姻对工资的影响从“0.1147”明显地下降到“0.0476”。尽管如此，这一结果在 99% 的水平上仍然是显著的，因此男人仍然可以预期婚姻会让收入增加。\n",
    "\n",
    "## 当面板数据对您没有帮助时\n",
    "\n",
    "使用面板数据和固定效应模型是一种非常强大的因果推理工具。当您没有随机数据也没有好的工具变量时，固定效应与非实验数据的因果推断一样令人信服。不过，值得一提的是，它并不是灵丹妙药。在某些情况下，即使面板数据也无济于事。\n",
    "\n",
    "最明显的一个情况是当您有随时间变化的混淆因素时。固定效应只能消除由于每个实体不变的属性带来的偏差。例如，假设您可以通过阅读书籍和食用大量优质脂肪来提高您的智力水平。同时这又会让你得到一份更高薪的工作和一个妻子。由于智力这个无法测量的混淆因素，固定效应将无法消除其带来的偏差，因为在这个例子中，智力是随时间变化的。\n",
    "\n",
    "![img](./data/img/fixed-effects/time-travel.png)\n",
    "\n",
    "固定效应失败的另一个不太明显的情况是当你有 **逆向因果（reversed causality）** 时。例如，假设不是婚姻让你赚得更多，而是当你收入越多，你结婚的机会就越多。在这种情况下，它们似乎具有正相关性，但是收入在线，是原因。二者都会随着时间在同一个方向发生变化，因此固定效应无法控制这一点。\n",
    "\n",
    "\n",
    "## 关键思想\n",
    "\n",
    "在这里，我们看到了如何使用面板数据，即我们在多个时间段内对同一个人进行多次测量的数据。在这种情况下，我们可以使用控制实体间差异的固定效应模型，保持所有个体不随时间变化的属性固定。这是一种强大且非常令人信服的控制混淆的方法，它可以得到与非随机数据一样好的结果。\n",
    "\n",
    "最后，我们看到 FE 并不是万能的。我们看到了两种不起作用的情况：当我们有反向因果关系时，以及当不可测量的混淆因素随时间变化时。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
